WGCNA Analysis

This is a Rmd file analyzing our raw count data by the WGCNA package as described in the manual and tutorials.

sessionInfo() #provides list of loaded packages and version of R. 
## R version 4.3.2 (2023-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.0
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.35     R6_2.5.1          fastmap_1.2.0     xfun_0.43        
##  [5] cachem_1.0.8      knitr_1.45        htmltools_0.5.8.1 rmarkdown_2.26   
##  [9] lifecycle_1.0.4   cli_3.6.2         sass_0.4.9        jquerylib_0.1.4  
## [13] compiler_4.3.2    rstudioapi_0.16.0 tools_4.3.2       evaluate_0.23    
## [17] bslib_0.7.0       yaml_2.3.8        rlang_1.1.3       jsonlite_1.8.8

First, download the necessary packages.

install.packages("BiocManager")
library("BiocManager")
BiocManager::install("impute", type = "source")
BiocManager::install("WGCNA",force = TRUE)
BiocManager::install("vsn")
install.packages("dendextend")

Next, load the packages

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(ggplot2)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(mgcv)
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
library(MuMIn)
library(emmeans)
library(locfit)
## locfit 1.5-9.9    2024-03-01
library(DESeq2)
## Loading required package: S4Vectors
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:lubridate':
## 
##     second, second<-
## The following object is masked from 'package:plyr':
## 
##     rename
## The following object is masked from 'package:tidyr':
## 
##     expand
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## The following object is masked from 'package:utils':
## 
##     findMatches
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following object is masked from 'package:lubridate':
## 
##     %within%
## The following object is masked from 'package:plyr':
## 
##     desc
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following object is masked from 'package:plyr':
## 
##     count
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
library(impute)
library(WGCNA)
## Loading required package: dynamicTreeCut
## Loading required package: fastcluster
## 
## Attaching package: 'fastcluster'
## The following object is masked from 'package:stats':
## 
##     hclust
## 
## 
## Attaching package: 'WGCNA'
## The following object is masked from 'package:IRanges':
## 
##     cor
## The following object is masked from 'package:S4Vectors':
## 
##     cor
## The following object is masked from 'package:stats':
## 
##     cor
library(xfun)
## 
## Attaching package: 'xfun'
## The following objects are masked from 'package:base':
## 
##     attr, isFALSE
library(genefilter)
## 
## Attaching package: 'genefilter'
## The following objects are masked from 'package:MatrixGenerics':
## 
##     rowSds, rowVars
## The following objects are masked from 'package:matrixStats':
## 
##     rowSds, rowVars
## The following object is masked from 'package:MASS':
## 
##     area
## The following object is masked from 'package:car':
## 
##     Anova
library(vsn)
library(RColorBrewer)
library(pheatmap)
library(dendextend)
## 
## ---------------------
## Welcome to dendextend version 1.17.1
## Type citation('dendextend') for how to cite the package.
## 
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
## 
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags: 
##   https://stackoverflow.com/questions/tagged/dendextend
## 
##  To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
## ---------------------
## 
## Attaching package: 'dendextend'
## The following object is masked from 'package:stats':
## 
##     cutree
sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.0
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] dendextend_1.17.1           pheatmap_1.0.12            
##  [3] RColorBrewer_1.1-3          vsn_3.68.0                 
##  [5] genefilter_1.82.1           xfun_0.43                  
##  [7] WGCNA_1.72-5                fastcluster_1.2.6          
##  [9] dynamicTreeCut_1.63-1       impute_1.74.1              
## [11] DESeq2_1.40.2               SummarizedExperiment_1.30.2
## [13] Biobase_2.60.0              MatrixGenerics_1.12.3      
## [15] matrixStats_1.3.0           GenomicRanges_1.52.1       
## [17] GenomeInfoDb_1.36.4         IRanges_2.34.1             
## [19] S4Vectors_0.38.2            BiocGenerics_0.46.0        
## [21] locfit_1.5-9.9              emmeans_1.10.1             
## [23] MuMIn_1.47.5                mgcv_1.9-1                 
## [25] nlme_3.1-164                MASS_7.3-60.0.1            
## [27] car_3.1-2                   carData_3.0-5              
## [29] ggplot2_3.5.1               lubridate_1.9.3            
## [31] Rmisc_1.5.1                 plyr_1.8.9                 
## [33] lattice_0.22-6              tidyr_1.3.1                
## [35] dplyr_1.1.4                
## 
## loaded via a namespace (and not attached):
##   [1] rstudioapi_0.16.0       jsonlite_1.8.8          magrittr_2.0.3         
##   [4] TH.data_1.1-2           estimability_1.5.1      rmarkdown_2.26         
##   [7] zlibbioc_1.46.0         vctrs_0.6.5             memoise_2.0.1          
##  [10] RCurl_1.98-1.14         base64enc_0.1-3         htmltools_0.5.8.1      
##  [13] S4Arrays_1.0.6          Formula_1.2-5           sass_0.4.9             
##  [16] bslib_0.7.0             htmlwidgets_1.6.4       sandwich_3.1-0         
##  [19] zoo_1.8-12              cachem_1.0.8            lifecycle_1.0.4        
##  [22] iterators_1.0.14        pkgconfig_2.0.3         Matrix_1.6-5           
##  [25] R6_2.5.1                fastmap_1.2.0           GenomeInfoDbData_1.2.10
##  [28] digest_0.6.35           colorspace_2.1-0        AnnotationDbi_1.62.2   
##  [31] Hmisc_5.1-2             RSQLite_2.3.6           fansi_1.0.6            
##  [34] timechange_0.3.0        httr_1.4.7              abind_1.4-5            
##  [37] compiler_4.3.2          bit64_4.0.5             withr_3.0.0            
##  [40] doParallel_1.0.17       htmlTable_2.4.2         backports_1.5.0        
##  [43] BiocParallel_1.34.2     viridis_0.6.5           DBI_1.2.2              
##  [46] DelayedArray_0.26.7     tools_4.3.2             foreign_0.8-86         
##  [49] nnet_7.3-19             glue_1.7.0              grid_4.3.2             
##  [52] checkmate_2.3.1         cluster_2.1.6           generics_0.1.3         
##  [55] gtable_0.3.5            preprocessCore_1.62.1   data.table_1.15.4      
##  [58] utf8_1.2.4              XVector_0.40.0          foreach_1.5.2          
##  [61] pillar_1.9.0            stringr_1.5.1           limma_3.56.2           
##  [64] splines_4.3.2           survival_3.6-4          bit_4.0.5              
##  [67] annotate_1.78.0         tidyselect_1.2.1        GO.db_3.17.0           
##  [70] Biostrings_2.68.1       knitr_1.45              gridExtra_2.3          
##  [73] stringi_1.8.3           yaml_2.3.8              evaluate_0.23          
##  [76] codetools_0.2-20        tibble_3.2.1            BiocManager_1.30.23    
##  [79] cli_3.6.2               affyio_1.70.0           rpart_4.1.23           
##  [82] xtable_1.8-4            munsell_0.5.1           jquerylib_0.1.4        
##  [85] Rcpp_1.0.12             coda_0.19-4.1           png_0.1-8              
##  [88] XML_3.99-0.16.1         parallel_4.3.2          blob_1.2.4             
##  [91] bitops_1.0-7            viridisLite_0.4.2       mvtnorm_1.2-5          
##  [94] scales_1.3.0            affy_1.78.2             purrr_1.0.2            
##  [97] crayon_1.5.2            rlang_1.1.3             KEGGREST_1.40.1        
## [100] multcomp_1.4-25

Load in raw count data

cts_raw <- read.csv("../../TagSeq_Output/HeronPdam_gene_count_matrix.csv") #load in data
rownames(cts_raw) <- cts_raw[,1] #set first column that contains gene names as rownames
cts_raw <- cts_raw[,-1] #remove the column with gene names
head(cts_raw)
##                                           RF13B_S85_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     3
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      2
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RF13D_S112_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      8
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RF14B_S116_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      4
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RF14C_S95_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          1
##                                           RF15B_S122_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      2
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RF15D_S84_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RF16A_S118_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RF16C_S97_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RF17B_S80_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                    10
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          1
##                                           RF17D_S99_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          1
##                                           RF18B_S90_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RF18D_S113_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RF19B_S103_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           1
##                                           RF19C_S83_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RF20B_S88_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     2
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RF20C_S115_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           1
##                                           RF22B_S91_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          1
##                                           RF22C_S126_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           1
##                                           RF23A_S125_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      5
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       6
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RF23C_S108_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RF24B_S92_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RF24D_S111_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       7
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RF25A_S96_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          1
##                                           RF25C_S110_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           1
##                                           RS11B_S86_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     1
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      6
## Pocillopora_acuta_HIv2___TS.g3788.t1                          1
##                                           RS11D_S87_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS12A_S105_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           1
##                                           RS12C_S106_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      2
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RS13A_S121_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RS13C_S119_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      2
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           1
##                                           RS14B_S107_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RS14C_S109_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      8
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RS15B_S104_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                      8
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                       0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                       0
## Pocillopora_acuta_HIv2___TS.g3788.t1                           0
##                                           RS15D_S82_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                    10
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS1B_S102_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     1
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      1
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS1C_S127_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     1
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     4
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          1
##                                           RS2B_S123_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS2C_S89_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                    0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                    0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                     0
## Pocillopora_acuta_HIv2___TS.g3788.t1                         1
##                                           RS3B_S100_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS3D_S117_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS6A_S81_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                    0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                    0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                     0
## Pocillopora_acuta_HIv2___TS.g3788.t1                         0
##                                           RS6D_S124_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     2
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS7B_S120_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS7C_S101_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS8B_S94_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                    0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                    5
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                     0
## Pocillopora_acuta_HIv2___TS.g3788.t1                         0
##                                           RS8C_S114_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                      0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                      0
## Pocillopora_acuta_HIv2___TS.g3788.t1                          0
##                                           RS9A_S98_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                    0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                    0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                    10
## Pocillopora_acuta_HIv2___TS.g3788.t1                         1
##                                           RS9C_S93_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1                    2
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1                    8
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1                     0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1                     5
## Pocillopora_acuta_HIv2___TS.g3788.t1                         1

Clean up sample names from “RF13B_S85_ALL.bam.gtf” to “RF13B”

colnames(cts_raw) <- gsub("_[^_]+$", "",colnames(cts_raw)) #get rid of "_ALL.bam.gtf"
colnames(cts_raw) <- gsub("_[^_]+$", "",colnames(cts_raw)) #get rid of "_S85"
head(colnames(cts_raw)) #see first 6 clean sample names
## [1] "RF13B" "RF13D" "RF14B" "RF14C" "RF15B" "RF15D"

Metadata from this dataset

coldata <- read.csv("../../../TagSeq_Submission/RNA Submission Sample List metadata.csv") #read in metadata file
coldata <- plyr::rename(coldata, c("Sample.Name"="Coral_ID")) #Make a column that represents the colonies
coldata$Colony <- gsub("A", "", coldata$Coral_ID) #ID which sample is from which colony by removing letters after colony name
coldata$Colony <- gsub("B", "", coldata$Colony) #ID which sample is from which colony by removing letters after colony name
coldata$Colony <- gsub("C", "", coldata$Colony) #ID which sample is from which colony by removing letters after colony name
coldata$Colony <- gsub("D", "", coldata$Colony) #ID which sample is from which colony by removing letters after colony name

coldata$Origin <- factor(coldata$Origin, levels = c("Slope","Flat")) #set variables to factors, with Slope as the baseline
coldata$Colony <- factor(coldata$Colony) #set variables to factors
coldata$Treatment <- factor(coldata$Treatment) #set variables to factors
head(coldata)
##   Coral_ID Origin Treatment Colony
## 1    RF13B   Flat  Variable   RF13
## 2    RF13D   Flat    Stable   RF13
## 3    RF14B   Flat  Variable   RF14
## 4    RF14C   Flat    Stable   RF14
## 5    RF15B   Flat  Variable   RF15
## 6    RF15D   Flat    Stable   RF15

Data sanity checks:

all(rownames(coldata$Coral_ID) %in% colnames(cts_raw)) #are all of the sample names (rows) in the metadata column names in the gene count matrix?
## [1] TRUE
all(rownames(coldata$Coral_ID) == colnames(cts_raw)) #are they the same in the same order?
## [1] TRUE

Filtering

dim(cts_raw)
## [1] 33730    48
cts_filt <- cts_raw[rowSums(!as.matrix(cts_raw)) < ncol(cts_raw), ] #here we remove all genes that were not expressed in any of our samples- we are left with 24,233 genes.

dim(cts_filt)
## [1] 24233    48

pOverA filtering to reduce dataset

ffun <- filterfun(pOverA(0.25,10))  #set up filtering parameters
cts_filtpoa <- genefilter((cts_filt), ffun) #apply filter
sum(cts_filtpoa) #count number of genes left
## [1] 9056
cts_filtpoa <- cts_filt[cts_filtpoa,] #keep only rows that passed filter, 9056 genes
dim(cts_filtpoa)
## [1] 9056   48
library("DESeq2")
dds <- DESeqDataSetFromMatrix(countData = cts_filtpoa,
                              colData=coldata,
                              design= ~Origin+Treatment)
dds
## class: DESeqDataSet 
## dim: 9056 48 
## metadata(1): version
## assays(1): counts
## rownames(9056): Pocillopora_acuta_HIv2___RNAseq.g27841.t1
##   Pocillopora_acuta_HIv2___RNAseq.g14011.t1 ...
##   Pocillopora_acuta_HIv2___RNAseq.g27413.t1
##   Pocillopora_acuta_HIv2___RNAseq.g10545.t1
## rowData names(0):
## colnames(48): RF13B RF13D ... RS9A RS9C
## colData names(4): Coral_ID Origin Treatment Colony
# transform the data 
vst <- vst(dds, blind=FALSE) # transform it vst

Effects of transformations on the variance

meanSdPlot(assay(vst))

Data quality assessment by sample clustering and visualization

Heatmap of the sample-to-sample distances

sampleDists <- dist(t(assay(vst)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vst$Origin, vst$Treatment, vst$Coral_ID, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors)

Principal component plot of the samples

withoutliers <- plotPCA(vst, intgroup=c("Origin", "Treatment"))
withoutliers

Removing outlier “RF16A” and “RF16C”

cts_raw_outrm <- cts_raw %>% dplyr::select(-c("RF16A","RF16C")) #remove those columns from count matrix
coldata_outrm <- coldata %>% filter(Coral_ID != c("RF16A","RF16C")) #removed those rows from metadata

Using filtered data instead of raw counts as input

cts_filt_outrm <- cts_raw_outrm[rowSums(!as.matrix(cts_raw_outrm)) < ncol(cts_raw_outrm), ] #here we remove all genes that were not expressed in any of our samples- we are left with 24,220 genes.

dim(cts_filt_outrm)
## [1] 24220    46

pOverA filtering to reduce dataset

ffun <- filterfun(pOverA(0.25,10))  #set up filtering parameters
cts_filt_outrm_poa <- genefilter((cts_filt_outrm), ffun) #apply filter
sum(cts_filt_outrm_poa) #count number of genes left
## [1] 9012
cts_filt_outrm_poa <- cts_filt_outrm[cts_filt_outrm_poa,] #keep only rows that passed filter, 9011 genes

cts_filtpoa <- cts_filt_outrm_poa

Data sanity checks:

all(rownames(coldata_outrm$Coral_ID) %in% colnames(cts_filt_outrm_poa)) #are all of the sample names (rows) in the metadata column names in the gene count matrix?
## [1] TRUE
all(rownames(coldata_outrm$Coral_ID) == colnames(cts_filt_outrm_poa)) #are they the same in the same order?
## [1] TRUE
write.csv(cts_filt_outrm_poa, "../../output/Filtered_gene_count_matrix.csv")
library("DESeq2")
dds2 <- DESeqDataSetFromMatrix(countData = cts_filtpoa,
                              colData=coldata_outrm,
                              design= ~Origin+Treatment)
dds2
## class: DESeqDataSet 
## dim: 9012 46 
## metadata(1): version
## assays(1): counts
## rownames(9012): Pocillopora_acuta_HIv2___RNAseq.g27841.t1
##   Pocillopora_acuta_HIv2___RNAseq.g14011.t1 ...
##   Pocillopora_acuta_HIv2___RNAseq.g27413.t1
##   Pocillopora_acuta_HIv2___RNAseq.g10545.t1
## rowData names(0):
## colnames(46): RF13B RF13D ... RS9A RS9C
## colData names(4): Coral_ID Origin Treatment Colony
# transform the data 
vst2 <- vst(dds2, blind=FALSE) # transform it vst

Effects of transformations on the variance

meanSdPlot(assay(vst2))

Data quality assessment by sample clustering and visualization

Heatmap of the sample-to-sample distances

sampleDists <- dist(t(assay(vst2)))
library("RColorBrewer")
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vst2$Origin, vst2$Treatment, vst2$Coral_ID, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors)

###Principal component plot of the samples

withoutoutliers <- plotPCA(vst2, intgroup=c("Origin", "Treatment"))
withoutoutliers

compare_figs <- cowplot::plot_grid(withoutliers, withoutoutliers,ncol=2, align="h")
compare_figs

Transpose the filtered gene count matrix so that the gene IDs are rows and the sample IDs are columns.

# transform the data 
vst2 <- vst(dds2) # transform it vst
vst2 <- assay(vst2) # call only the transformed coutns in the dds object
vst2 <- t(vst2) # transpose columns to rows and vice versa

Check for genes and samples with too many missing values with goodSamplesGenes. There shouldn’t be any because we performed pre-filtering

dim(vst2) #  9012 genes; 46  samples
## [1]   46 9012
gsg <- goodSamplesGenes(vst2, verbose = 3);  # We first check for genes and samples with too many missing values
##  Flagging genes and samples with too many missing values...
##   ..step 1
gsg$allOK # If the statement returns TRUE, all genes have passed the cuts. 
## [1] TRUE

Conduct WGCNA

library(WGCNA)
# The following setting is important, do not omit.
options(stringsAsFactors = FALSE)
##Soft threshold
dim(vst2) #  46 9012
## [1]   46 9012
# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function

#the below takes a long time to run, so is commented out and the pre-run results are loaded in below
## sft = pickSoftThreshold(vst2, powerVector = powers, verbose = 5) #...wait for this to finish
## save(sft, file = "../../output/WGCNA/sft.RData")
load("../../output/WGCNA/sft.RData")

# pickSoftThreshold 
#  performs the analysis of network topology and aids the
# user in choosing a proper soft-thresholding power.
# The user chooses a set of candidate powers (the function provides suitable default values)
# function returns a set of network indices that should be inspected

sizeGrWindow(9, 5) # set window size 
# png to output 
png("../../output/WGCNA/sft.png", 1000, 1000, pointsize=20)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));

text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");

# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red") 

# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
dev.off() # output 
## pdf 
##   2
#I used a scale-free topology fit index **R^2 of 0.9**. This lowest recommended R^2 by Langfelder and Horvath is 0.8. I chose 0.9 because we want to use the smallest soft thresholding power that maximizes with model fit. It appears that our **soft thresholding power is 5** because it is the lowest power above the R^2=0.9 threshold that maximizes with model fit.  

Look for outliers by examining tree of samples

sampleTree = hclust(dist(vst2), method = "average") # Next we cluster the samples (in contrast to clustering genes that will come later)  to see if there are any obvious outliers.There don't look to be any outliers, so we will move on with business as usual.  
sizeGrWindow(12,9) 
par(cex = 0.6);
par(mar = c(0,4,2,0))
plot(sampleTree)

Start the step-wise module construction:

Step 1: Create adjacency matrix

softPower = 5 # set the soft threshold based on the plots above 

# signed 
#to get this to work with a ton of genes >30K, I had to increase my memory limit using: memory.limit(size = 35000) 
adjacency_sign = adjacency(vst2, power = softPower, type="signed")  #Calculate adjacency

Step 2: Turn adjacency into topological overlap: Calculation of the topological overlap matrix, (TOM) and the corresponding dissimilarity, from a given adjacency matrix.

#the below takes a long time to run, so is commented out and the pre-run results are loaded in below
## TOM_sign = TOMsimilarity(adjacency_sign, TOMType="signed") #Translate adjacency into topological overlap matrix
## save(TOM_sign, file = "../../output/WGCNA/TOM_sign.Rdata")
load("../../output/WGCNA/TOM_sign.Rdata")

dissTOM_sign   = 1-TOM_sign

Step 3: Call the hierarchical clustering function - plot the tree

# Call the hierarchical clustering function
#to get this to work, I had to increase my memory limit using: memory.limit(size = 45000) 
geneTree_sign   = hclust(as.dist(dissTOM_sign), method = "average");

# Plot the resulting clustering tree (dendrogram) Each leaf corresponds to a gene, branches grouping together densely are interconnected, highly co-expressed genes.  
sizeGrWindow(12,9)

plot(geneTree_sign, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity - SIGNED",
     labels = FALSE, hang = 0.04);

Step 4: Set module size and ‘cutreeDynamic’ to create clusters

#Module identification is essentially cutting the branches off the tree in the dendrogram above. We like large modules, so we set the **minimum module size** relatively high, so we will set the minimum size at 30. 
minModuleSize = 30; # set this for the subseqent call...

dynamicMods_sign = cutreeDynamic(dendro = geneTree_sign, distM = dissTOM_sign,
                            deepSplit = 1, pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize);
##  ..cutHeight not given, setting it to 0.959  ===>  99% of the (truncated) height range in dendro.
##  ..done.
table(dynamicMods_sign) # number of genes per module. Module 0 is reserved for unassigned genes. The are other modules will be listed largest to smallest. 
## dynamicMods_sign
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
## 2558 1809  942  917  500  425  396  220  219  218  180  170  154  105   71   65 
##   17 
##   63

Step 5: convert numeric network to colors and plot the dendrogram

# Convert numeric lables into colors
dynamicColors_sign = labels2colors(dynamicMods_sign) # add colors to module labels (previously numbers)
table(dynamicColors_sign) # lets look at this table...
## dynamicColors_sign
##        black         blue        brown         cyan        green  greenyellow 
##          396         1809          942          105          500          180 
##       grey60    lightcyan      magenta midnightblue         pink       purple 
##           63           65          219           71          220          218 
##          red       salmon          tan    turquoise       yellow 
##          425          154          170         2558          917
# Plot the dendrogram and colors underneath

plotDendroAndColors(geneTree_sign, dynamicColors_sign, "Dynamic Tree Cut - SIGNED",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05,
                    main = "Gene dendrogram and module colors 'SIGNED'")

sizeGrWindow(8,6)
png("../../output/WGCNA/GeneDendrogram.png", 1000, 1000, pointsize=20)
plotDendroAndColors(geneTree_sign, dynamicColors_sign, "Dynamic Tree Cut - SIGNED",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05,
                    main = "Gene dendrogram and module colors 'SIGNED'")
dev.off()
## pdf 
##   2

Step 6: Calculate Eigengenes - view thier connectivity based on ‘MEDiss = 1-cor(MEs)’

# Calculate eigengenes
MEList = moduleEigengenes(vst2, colors = dynamicColors_sign, softPower = 5)

MEs = MEList$eigengenes
MEs
##            MEblack        MEblue      MEbrown        MEcyan     MEgreen
## RF13B  0.035738890 -0.0005614067  0.086214000  0.0230457125  0.06457070
## RF13D  0.016715232  0.0355808790  0.077795819  0.0108599933  0.04977675
## RF14B -0.146566530 -0.0989830172  0.243945678  0.0089738153  0.09690707
## RF14C -0.144417844  0.0293510325  0.231391380  0.0555552483  0.12492062
## RF15B  0.361363331  0.0125541542  0.195194429  0.0669848895  0.09848360
## RF15D  0.340691436 -0.0628815348  0.130600265  0.0613695551  0.15041068
## RF17B -0.166389167 -0.0083085168  0.198989474 -0.0330477836  0.06580995
## RF17D -0.174336507  0.0274499964  0.242761454  0.0109581641  0.20156729
## RF18B -0.163062969  0.0730460426  0.213165173  0.0428568437  0.15059016
## RF18D -0.198796035  0.1501095732  0.221598881  0.0156327910  0.26172808
## RF19B -0.052085383  0.1809702199  0.156245872  0.0534562058  0.11871425
## RF19C -0.133192276 -0.1908284494  0.207082710  0.0058409291  0.09057738
## RF20B -0.203191630  0.1418285398  0.167736368  0.0175489549  0.09355805
## RF20C -0.079445376  0.0589582967  0.220340540  0.0517627628  0.08359879
## RF22B  0.013623481 -0.1171079552  0.003985773  0.0146616850 -0.06040064
## RF22C -0.011991191  0.1497873255  0.042606148  0.0054811614 -0.06439484
## RF23A  0.312454266  0.1092940451  0.061065835  0.0347675690  0.09520240
## RF23C  0.295349869  0.0345798457  0.114909745  0.0202737297  0.03245128
## RF24B  0.356715054  0.1929708845  0.106044589  0.0913070182  0.03457393
## RF24D  0.306302794  0.1516425029  0.076269408  0.0729544308  0.08671047
## RF25A  0.080575924 -0.0834460336 -0.007865851 -0.0080197431 -0.04258953
## RF25C  0.093458065  0.1069453692 -0.019058548  0.0240384357  0.09631201
## RS11B -0.004644330 -0.1099408605 -0.092245924  0.0376314369  0.10628308
## RS11D -0.026519519  0.0766325673 -0.153945434  0.0009392144  0.34571416
## RS12A -0.018592944  0.0934519345 -0.215804220  0.0079270236 -0.25126870
## RS12C  0.025448390 -0.2497580080 -0.180395235  0.0528064649 -0.06947225
## RS13A -0.013832874  0.2242246956 -0.087780763  0.0043265989 -0.17859214
## RS13C -0.016960538  0.0901775914 -0.092445893  0.0264009339 -0.20267882
## RS14B -0.025815615 -0.1101678607 -0.138277199  0.0152473895 -0.14673511
## RS14C -0.093490223 -0.3533046553 -0.224348383 -0.0853239910 -0.29804769
## RS15B -0.043303446 -0.0724784483 -0.135120081  0.0254364294 -0.25047579
## RS15D -0.063247364 -0.4464606863 -0.250028576  0.0434992663 -0.21085363
## RS1B  -0.071734949  0.1666814228 -0.146538808  0.0318043922  0.01154907
## RS1C  -0.081189863 -0.0208045578 -0.082597955  0.0382703168  0.03230402
## RS2B   0.081141470  0.0226962013 -0.046649316 -0.0191289688 -0.10542840
## RS2C   0.001196774  0.0341686249 -0.061856762  0.0202354959 -0.13064096
## RS3B  -0.007344414  0.0096703526 -0.121890691  0.0118973832 -0.01743041
## RS3D  -0.005781953  0.0046080335 -0.079855011  0.0442362556  0.03886139
## RS6A  -0.070105162 -0.0231943131 -0.147206872 -0.9641101677 -0.14047140
## RS6D  -0.055100416 -0.1924917784 -0.048878771 -0.0740728504 -0.01160300
## RS7B  -0.115318529 -0.0362586532 -0.189674392 -0.0186040960 -0.26884057
## RS7C  -0.191763564 -0.3670929579 -0.153268154  0.0157100485 -0.22995477
## RS8B   0.035698058  0.1335243682 -0.069128256  0.0594849913 -0.01327168
## RS8C   0.017104694  0.1220210827 -0.089764315  0.0442748956 -0.03140854
## RS9A  -0.008609524  0.0942830753 -0.113018368  0.0003933584  0.23592747
## RS9C   0.013252409  0.0168610358 -0.050299764  0.0334558101 -0.04254375
##       MEgreenyellow     MEgrey60  MElightcyan    MEmagenta MEmidnightblue
## RF13B    0.03096313  0.048813655  0.004539119  0.040533178  -0.0245491734
## RF13D    0.07466948 -0.002747228 -0.048921622  0.004749849  -0.0183763055
## RF14B   -0.02351578  0.052054548  0.007711938 -0.055936540  -0.0228806751
## RF14C    0.05431928  0.103521196 -0.069472828 -0.015987095   0.0063273755
## RF15B    0.11598069 -0.027006621 -0.344930466 -0.362117043   0.1513293329
## RF15D   -0.03629918  0.005716079 -0.167570877 -0.295093407   0.2593929751
## RF17B    0.04043036  0.037559058 -0.061350955 -0.105576226   0.0522009454
## RF17D    0.08347539  0.240079479 -0.003870975 -0.071316708   0.0889135020
## RF18B    0.17271940  0.088913590 -0.120851396 -0.157907924   0.0770626674
## RF18D    0.16454055  0.312762543 -0.045673177 -0.139895717   0.1822861472
## RF19B    0.11959708  0.096034394  0.029738991 -0.021384698  -0.0659116912
## RF19C   -0.06142188  0.035342775 -0.226848538 -0.113708482   0.0063878362
## RF20B    0.20131380  0.085070425 -0.066774285 -0.071679981  -0.0725077225
## RF20C    0.03843045  0.065751826 -0.005878656 -0.079070138  -0.0267314858
## RF22B   -0.12513476 -0.112608524  0.034288078 -0.033064166  -0.0109426095
## RF22C    0.18911651 -0.134875929  0.010881417 -0.037849207  -0.1930922669
## RF23A    0.07443417 -0.015840141 -0.080743552 -0.255329622   0.1073787126
## RF23C    0.12020491 -0.135276896 -0.256935676 -0.350307537   0.0506762137
## RF24B    0.16937663 -0.090056137 -0.119388795 -0.241229366   0.0014127112
## RF24D    0.11689102  0.002523370  0.002249586 -0.107578383   0.0916879301
## RF25A   -0.11865150 -0.047263339  0.136314599  0.060319789   0.0447707556
## RF25C    0.06960809  0.065371487  0.115822838  0.091429118   0.1221668700
## RS11B   -0.12420996  0.061302028 -0.150987001  0.022147377   0.1402703642
## RS11D    0.05539636  0.541732945 -0.060301638 -0.009668998   0.4557546749
## RS12A   -0.07207071 -0.112117851  0.410469539  0.213579019  -0.2390647714
## RS12C   -0.29590257 -0.015964786  0.110170934  0.097511706   0.0886475060
## RS13A    0.13529317 -0.157437914  0.243742958  0.153568478  -0.3479725582
## RS13C    0.05040341 -0.225760997  0.136072703  0.120176386  -0.2901638771
## RS14B   -0.17606312 -0.134924778  0.149095999  0.126931780  -0.1067129932
## RS14C   -0.32994363 -0.250860431  0.111490397  0.110626775  -0.1183388272
## RS15B   -0.10615944 -0.212941292  0.118437262  0.196764982  -0.2220598840
## RS15D   -0.50142342 -0.122195112 -0.002867751  0.192810781   0.0205032779
## RS1B     0.09307705  0.102779430  0.092266234  0.165333200  -0.0786940195
## RS1C     0.03100818  0.073828989 -0.122813217  0.021896321   0.1067984725
## RS2B     0.03704151 -0.153620437 -0.026877271 -0.043371606  -0.0752202101
## RS2C     0.02758929 -0.078462245  0.011979459  0.049463264   0.0101279467
## RS3B     0.02123936 -0.033816965 -0.028347890 -0.084532348  -0.0601723773
## RS3D     0.01938915  0.033029583 -0.081046733 -0.004556133  -0.0001098878
## RS6A     0.07246209 -0.073104271 -0.029397698  0.039469343  -0.1116858276
## RS6D    -0.04656323  0.071344990 -0.282674982  0.034322064   0.1192637551
## RS7B    -0.13664499 -0.108397584  0.272633413  0.224680179  -0.1356225698
## RS7C    -0.36398426 -0.200868078 -0.047441111  0.181188310   0.0559780216
## RS8B     0.03695795  0.064064360  0.197500072  0.179495826  -0.1395988621
## RS8C     0.03472902  0.045298235  0.277202598  0.222178199  -0.1194076648
## RS9A     0.08963872  0.288165080 -0.078615357 -0.010370655   0.2943785816
## RS9C    -0.02230774 -0.074912507  0.057974312  0.118356057  -0.0539003158
##              MEpink      MEpurple         MEred     MEsalmon        MEtan
## RF13B -0.0647255585  0.0151669067 -0.0412048930 -0.029662215  0.004223004
## RF13D  0.0567711235  0.0045436211  0.0182869908  0.049450663 -0.066235131
## RF14B -0.0227504438  0.2937518176  0.0275893621  0.103027496  0.244085251
## RF14C -0.0264419838  0.2502514788  0.0597864506  0.001006515  0.166964022
## RF15B  0.2977573509 -0.1588172587  0.2684242914  0.307372343 -0.248831159
## RF15D  0.2675751469 -0.1318729337  0.2400972161  0.345022730 -0.142520988
## RF17B  0.0201348420  0.2433856181  0.0331804376  0.081061125  0.135198540
## RF17D -0.0450480509  0.3059433164  0.0755925194 -0.096919085  0.250243406
## RF18B  0.1219523070  0.2395830981  0.1781091787 -0.020943218  0.091557210
## RF18D  0.0163742929  0.2744014056  0.0723252141 -0.127858562  0.059736901
## RF19B -0.1720667112  0.2498582943 -0.0512888820 -0.108542763  0.116936546
## RF19C  0.1635415194  0.2285495607  0.1172686495  0.229733153  0.163346305
## RF20B  0.0648717076  0.2523383762  0.1687415292 -0.114089695  0.103224721
## RF20C -0.0006323075  0.2815220169  0.1134350350  0.015289757  0.220902418
## RF22B  0.0149928171  0.0064553830  0.0784595688  0.066325514  0.018193497
## RF22C  0.0863902854  0.0148348145  0.2755034633 -0.231270487 -0.002092525
## RF23A  0.0596190520 -0.1409510451  0.1214340863  0.105376690 -0.210702052
## RF23C  0.3621427737 -0.1708847815  0.3540544053  0.288780993 -0.283396351
## RF24B  0.0647018282 -0.1184864108  0.0693294328  0.104284706 -0.251013442
## RF24D  0.0353817295 -0.1166670846  0.1366884407  0.031922256 -0.114720296
## RF25A -0.0562018151  0.0289125326 -0.0087538500  0.068680685  0.067021221
## RF25C -0.1662642378  0.0155097767 -0.1071282788 -0.140434773  0.027334367
## RS11B  0.1710590129 -0.0851167145  0.0177296440  0.071706946 -0.079916232
## RS11D -0.0123017550 -0.0628538011 -0.0062733166 -0.076054829 -0.015876831
## RS12A -0.4251690594 -0.0525968459 -0.2413145580 -0.199216173  0.121356930
## RS12C -0.2233364331 -0.0952775719 -0.1958107921  0.158852236  0.209851744
## RS13A -0.0548748042  0.0002636796 -0.1072797695 -0.183481687 -0.095731690
## RS13C  0.1051871651 -0.0270560269  0.0968598729 -0.099625983  0.009833324
## RS14B -0.0590788210 -0.1035488583 -0.1273157643  0.042699790  0.078095475
## RS14C  0.0352559251 -0.1717060997 -0.1259982955  0.097974513  0.114293125
## RS15B  0.0304411988 -0.0565690422 -0.0572671524 -0.067166091  0.116585393
## RS15D -0.2099670058 -0.0853186762 -0.3806235298  0.071515739  0.114460742
## RS1B  -0.2663953476 -0.1157331594 -0.2840088869 -0.330538257 -0.246867986
## RS1C   0.0357629036 -0.1631654026 -0.0767202606  0.029871371 -0.237046609
## RS2B   0.0606121877 -0.0915451370 -0.0867856299  0.140848204 -0.128993836
## RS2C   0.0099759699 -0.1103592577 -0.0290874453 -0.012361104 -0.183766930
## RS3B   0.0387810590 -0.0995603565  0.0381724975  0.007501047 -0.204687975
## RS3D   0.0481439887 -0.0967228437  0.0557197753 -0.014781091 -0.124814400
## RS6A   0.0257448243 -0.0793110778 -0.0310425919 -0.051437010 -0.088565740
## RS6D   0.1257730641 -0.0832343533 -0.0373659007  0.129638467 -0.088335417
## RS7B  -0.1065909184 -0.0385199251 -0.1463727209 -0.183028241  0.161101074
## RS7C   0.0259720871 -0.1106157481 -0.1274810919  0.158323826  0.151502090
## RS8B  -0.2639940985 -0.0451161373 -0.1818321106 -0.202533304  0.102274828
## RS8C  -0.2604526033 -0.0217115573 -0.1850355018 -0.267425678  0.045501922
## RS9A   0.0598996333 -0.0546996378  0.0188980123 -0.093487238 -0.062180990
## RS9C   0.0314761592 -0.0172539519  0.0003051485 -0.055409283 -0.017527481
##        MEturquoise     MEyellow
## RF13B -0.049987446  0.063928602
## RF13D -0.019137147 -0.027995593
## RF14B -0.093185651  0.087636966
## RF14C -0.147709679  0.028957143
## RF15B  0.035232756 -0.382862198
## RF15D  0.078786749 -0.140150323
## RF17B -0.103854472  0.071280185
## RF17D -0.215289310  0.133312986
## RF18B -0.192378713  0.007020303
## RF18D -0.265041143  0.169001768
## RF19B -0.186757003  0.128738799
## RF19C  0.040882662 -0.155862228
## RF20B -0.249292334  0.114648118
## RF20C -0.150262215  0.146705924
## RF22B  0.091541606  0.028452390
## RF22C -0.186848856  0.090375080
## RF23A -0.089565674  0.014122502
## RF23C  0.013769414 -0.153623198
## RF24B -0.105806483 -0.039796849
## RF24D -0.105876790 -0.006773177
## RF25A  0.040721698  0.132236412
## RF25C -0.081198520  0.122990690
## RS11B  0.130518859 -0.158996790
## RS11D -0.014869673  0.054816786
## RS12A -0.007885866  0.293024916
## RS12C  0.308619894 -0.086322820
## RS13A -0.147909312  0.216107093
## RS13C -0.027925323  0.115750928
## RS14B  0.165844629  0.057513137
## RS14C  0.326288209  0.001611669
## RS15B  0.130767735 -0.019966361
## RS15D  0.307123864 -0.172427392
## RS1B  -0.101737577  0.104667302
## RS1C   0.143128087 -0.183411878
## RS2B   0.100982719 -0.083769236
## RS2C   0.024999718 -0.048380219
## RS3B   0.024856931 -0.037476547
## RS3D   0.065598226 -0.130456228
## RS6A   0.076021058 -0.120889165
## RS6D   0.237997729 -0.426863340
## RS7B   0.010898755  0.148428399
## RS7C   0.300330368 -0.232336901
## RS8B  -0.059359541  0.115675577
## RS8C  -0.069784904  0.204699604
## RS9A  -0.029634459 -0.091637783
## RS9C   0.046386423  0.048294943
library(dplyr)
MEs <- MEs %>% 
    select_if(~ !any(is.na(.)))

# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs, use = "pairwise.complete.obs");

# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result
sizeGrWindow(7, 6)
png("../../output/WGCNA/ClusterEigengenes.png", 1000, 1000, pointsize=20)
plot(METree, main = "Clustering of module eigengenes - SIGNED (dissimilarity calc = MEDiss = 1-cor(MEs))",
     xlab = "", sub = "")
MEDissThres = 0.15 
abline(h=MEDissThres, col = "red")
dev.off()
## pdf 
##   2

Step 7: Specify the cut line for the dendrogram (module) - Calc MODULE EIGENGENES (mergeMEs)

We had 17 modules before merging, and 15 modules after merging

MEDissThres = 0.15 # **Merge modules with >85% eigengene similarity.** Most studies use somewhere between 80-90% similarity. I will use 85% similarity as my merging threshold.
# Plot the cut line into the dendrogram
#abline(h=MEDissThres, col = "red")
# Call an automatic merging function
# merge = mergeCloseModules(dds.d0_vst, dynamicColors, cutHeight = MEDissThres, verbose = 3)
merge = mergeCloseModules(vst2, dynamicColors_sign, cutHeight = MEDissThres, verbose = 3)
##  mergeCloseModules: Merging modules whose distance is less than 0.15
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 17 module eigengenes in given set.
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 15 module eigengenes in given set.
##    Calculating new MEs...
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 15 module eigengenes in given set.
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEs;
library(dplyr)
mergedMEs <- mergedMEs %>% 
    select_if(~ !any(is.na(.)))
# Cluster module eigengenes
MEDiss2 = 1-cor(mergedMEs,use = 'pairwise.complete.obs');
MEDiss2
##                   MEcyan   MEblack     MEsalmon    MEpink     MEred   MEgreen
## MEcyan         0.0000000 0.8363810 9.456342e-01 1.0357865 0.9335928 0.8183041
## MEblack        0.8363810 0.0000000 5.863287e-01 0.6853443 0.6615370 0.9676272
## MEsalmon       0.9456342 0.5863287 3.330669e-16 0.3469933 0.5950190 0.9701274
## MEpink         1.0357865 0.6853443 3.469933e-01 0.0000000 0.1847992 0.7966365
## MEred          0.9335928 0.6615370 5.950190e-01 0.1847992 0.0000000 0.5869122
## MEgreen        0.8183041 0.9676272 9.701274e-01 0.7966365 0.5869122 0.0000000
## MEmidnightblue 0.8698516 0.8018912 5.817853e-01 0.6693067 0.7309443 0.2592991
## MEblue         0.9280028 0.7794428 1.490247e+00 1.0050357 0.6357651 0.5714618
## MEbrown        0.8001074 0.9517543 7.364483e-01 0.5916456 0.3078619 0.4095681
## MEpurple       0.9256202 1.5772808 1.197132e+00 1.0905539 0.7866850 0.5726397
## MEturquoise    1.1349435 0.9426693 5.237657e-01 0.9541162 1.4074916 1.5287035
## MEtan          0.9540523 1.6478878 1.185038e+00 1.4159836 1.2577153 1.0651031
## MEmagenta      1.0943541 1.4972583 1.596656e+00 1.7187830 1.8343444 1.5388426
## MElightcyan    0.9911347 1.2290197 1.671253e+00 1.7962628 1.6277834 1.4917147
## MEyellow       0.8784473 1.2609555 1.706417e+00 1.5933035 1.2003612 1.0077827
##                MEmidnightblue    MEblue   MEbrown  MEpurple MEturquoise
## MEcyan              0.8698516 0.9280028 0.8001074 0.9256202   1.1349435
## MEblack             0.8018912 0.7794428 0.9517543 1.5772808   0.9426693
## MEsalmon            0.5817853 1.4902471 0.7364483 1.1971317   0.5237657
## MEpink              0.6693067 1.0050357 0.5916456 1.0905539   0.9541162
## MEred               0.7309443 0.6357651 0.3078619 0.7866850   1.4074916
## MEgreen             0.2592991 0.5714618 0.4095681 0.5726397   1.5287035
## MEmidnightblue      0.0000000 1.0911109 0.7673563 1.0220637   0.9399263
## MEblue              1.0911109 0.0000000 0.6094019 0.7720365   1.8569108
## MEbrown             0.7673563 0.6094019 0.0000000 0.2973886   1.6531319
## MEpurple            1.0220637 0.7720365 0.2973886 0.0000000   1.6328755
## MEturquoise         0.9399263 1.8569108 1.6531319 1.6328755   0.0000000
## MEtan               1.1727525 1.3358764 0.8763249 0.3437970   1.0458386
## MEmagenta           1.4884948 1.2900461 1.7130179 1.1285645   0.6777335
## MElightcyan         1.6092707 0.9317914 1.4922356 1.0135890   1.0509252
## MEyellow            1.4254228 0.5307596 0.9419223 0.5606714   1.6026729
##                    MEtan MEmagenta  MElightcyan  MEyellow
## MEcyan         0.9540523 1.0943541 9.911347e-01 0.8784473
## MEblack        1.6478878 1.4972583 1.229020e+00 1.2609555
## MEsalmon       1.1850376 1.5966556 1.671253e+00 1.7064168
## MEpink         1.4159836 1.7187830 1.796263e+00 1.5933035
## MEred          1.2577153 1.8343444 1.627783e+00 1.2003612
## MEgreen        1.0651031 1.5388426 1.491715e+00 1.0077827
## MEmidnightblue 1.1727525 1.4884948 1.609271e+00 1.4254228
## MEblue         1.3358764 1.2900461 9.317914e-01 0.5307596
## MEbrown        0.8763249 1.7130179 1.492236e+00 0.9419223
## MEpurple       0.3437970 1.1285645 1.013589e+00 0.5606714
## MEturquoise    1.0458386 0.6777335 1.050925e+00 1.6026729
## MEtan          0.0000000 0.6018516 6.027920e-01 0.5881998
## MEmagenta      0.6018516 0.0000000 2.388853e-01 0.6635948
## MElightcyan    0.6027920 0.2388853 1.110223e-16 0.2320545
## MEyellow       0.5881998 0.6635948 2.320545e-01 0.0000000
METree2 = hclust(as.dist(MEDiss2), method = "average");
# Plot the result
plot(METree2, main = "Clustering of module eigengenes - SIGNED (dissimilarity calc = MEDiss = 1-cor(MEs))",
     xlab = "", sub = "")

sizeGrWindow(7, 6)
png("../../output/WGCNA/ClusterEigengenes_merged.png", 1000, 1000, pointsize=20)
plot(METree2, main = "Clustering of module eigengenes - SIGNED (dissimilarity calc = MEDiss = 1-cor(MEs))",
     xlab = "", sub = "")
dev.off()
## pdf 
##   2

Step 8: Plot dendrogram with the cut line ‘MEDissThres’

plotDendroAndColors(geneTree_sign, cbind(dynamicColors_sign, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

sizeGrWindow(12, 9)

png("../../output/WGCNA/ClusterDendrogram_signed.png", 1000, 1000, pointsize=20)
plotDendroAndColors(geneTree_sign, cbind(dynamicColors_sign, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
dev.off()
## pdf 
##   2

Step 9: Commit to mergedcolors as ‘MEs’ and ‘moduleColors’

# Rename to moduleColors
moduleColors = mergedColors
# Construct numerical labels corresponding to the colors
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
MEs = mergedMEs;
# Save module colors and labels for use in subsequent parts
save(MEs, moduleLabels, moduleColors, geneTree_sign, file = "../../output/WGCNA/networkConstruction-stepByStep.RData")
# write csv - save the module eigengenes
write.csv(MEs, file = "../../output/WGCNA/WGCNA_ModuleEigengenes.csv")
table(mergedColors)
## mergedColors
##        black         blue        brown         cyan        green    lightcyan 
##          396         1989          942          105          563           65 
##      magenta midnightblue         pink       purple          red       salmon 
##          219           71          220          218          425          154 
##          tan    turquoise       yellow 
##          170         2558          917

Prepare for module trait associations - Eigengene calc - trait data as factors

#Prepare trait data. Data has to be numeric, so I will substitute time_points and type for numeric values. The "trait" we are considering here is habitat of origin. Make a dataframe that has a column for each origin name and a row for samples. Populate a 1 for samples that match each origin = flat and a 0 for samples not matching (e.g., origin = slope). This process changes origin from a categorical variable into a binary variable. This will allow for correlations between mean eigengenes and habitat of origin.  
Samples = rownames(vst2);# start new variable 'd21.Samples' calling the row names of the gene data (sample as 'Mouse' in trait data)
TreatRows = match(Samples, coldata$Coral_ID); # match the names 
Traits = coldata[TreatRows, -1]; # removes the row numbers
rownames(Traits) = coldata[TreatRows, 1]; # inserts the new traitRows - matches sample treatments
dim(Traits) #  46 Samples 2 columns (Origin and Treatment)
## [1] 46  3
all(rownames(Traits) == rownames(vst2))  # should be TRUE
## [1] TRUE
dim(Traits) #  46 2
## [1] 46  3

Origin

# ALL TRAITS 
head(Traits)
##       Origin Treatment Colony
## RF13B   Flat  Variable   RF13
## RF13D   Flat    Stable   RF13
## RF14B   Flat  Variable   RF14
## RF14C   Flat    Stable   RF14
## RF15B   Flat  Variable   RF15
## RF15D   Flat    Stable   RF15
# Origin groups 
Traits.Origin <-  Traits %>% dplyr::select('Origin')

Traits.Origin$Flat <- (Traits.Origin$Origin == "Flat")
Traits.Origin$Flat   <- as.numeric(Traits.Origin$Flat)
Traits.Origin$Flat <- as.factor(Traits.Origin$Flat)

Traits.Origin$Slope <- (Traits.Origin$Origin == "Slope")
Traits.Origin$Slope   <- as.numeric(Traits.Origin$Slope)
Traits.Origin$Slope <- as.factor(Traits.Origin$Slope)

Traits.Origin <- Traits.Origin[,c(2:3)]
head(Traits.Origin)
##       Flat Slope
## RF13B    1     0
## RF13D    1     0
## RF14B    1     0
## RF14C    1     0
## RF15B    1     0
## RF15D    1     0

Cluster samples by Origin

# Primary treatment Only
traitColors_Primary = labels2colors(Traits.Origin); # Convert traits to a color representation: white means low, red means high, grey means missing entry
plotDendroAndColors(sampleTree, traitColors_Primary, # Plot the sample dendrogram and the colors underneath.
                    groupLabels = names(Traits.Origin), 
                    main = "Sample dendrogram and trait heatmap")

png("../../output/WGCNA/ClusterTree_Origin.png", 1000, 1000, pointsize=20)
traitColors_Primary = labels2colors(Traits.Origin); # Convert traits to a color representation: white means low, red means high, grey means missing entry
plotDendroAndColors(sampleTree, traitColors_Primary, # Plot the sample dendrogram and the colors underneath.
                    groupLabels = names(Traits.Origin), 
                    main = "Sample dendrogram and trait heatmap")
dev.off()
## quartz_off_screen 
##                 2
# save data
save(vst2, Traits, file = "../../output/WGCNA/dataInput.RData")

# write the vst transformed data 
write.csv(vst2, "../../output/WGCNA/vstTransformed_WGCNAdata.csv") # write
# identify modules that are significantly associated with the measured clinical traits.

# Since we already have a summary profile (eigengene) for each module, we simply correlate eigengenes with external traits and look for the most significant associations:

# Define numbers of genes and samples
nGenes = ncol(vst2); # 9266
nSamples = nrow(vst2); # 46
# # Recalculate MEs with color labels
MEs0 = moduleEigengenes(vst2, moduleColors)$eigengenes
MEs = orderMEs(MEs0) # reorders the columns (colors/modules) , doesn't change anything

# change character treatments to integers
# ALL TRAIT DATA
Traits$Origin <- as.factor(Traits$Origin)
Traits$Origin <- as.numeric(Traits$Origin)

Module trait correlation

# ALL TRAIT DATA

dim(Traits)  # 48 2
## [1] 46  3
dim(MEs)  # 48 15
## [1] 46 15
# moduleTraitCor = cor(MEs, d0.Traits, use = "p");
# moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

nSamples = nrow(vst2)

# Origin 
Traits.Origin$Flat <- as.numeric(Traits.Origin$Flat)
Traits.Origin$Slope <- as.numeric(Traits.Origin$Slope)
moduleTraitCor_Primary = cor(MEs, Traits.Origin, use = "p");
moduleTraitPvalue_Primary = corPvalueStudent(moduleTraitCor_Primary, nSamples);

Heatmaps

sizeGrWindow(10,10)
# Will display correlations and their p-values
d0.PRIMARYTreatments.matrix <-  paste(signif(moduleTraitCor_Primary, 3), "\n(",
                                       signif(moduleTraitPvalue_Primary, 3), ")", sep = "")
#dim(textMatrix) == dim(moduleTraitCor_treatonly)
par(mar = c(8, 9.5, 5, 3));
# Display the correlation values within a heatmap plot

labeledHeatmap(Matrix = moduleTraitCor_Primary,
               xLabels = names(Traits.Origin),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = TRUE,
               colors = blueWhiteRed(50),
               textMatrix = d0.PRIMARYTreatments.matrix,
               setStdMargins = FALSE,
               cex.text = 1,
               zlim = c(-1,1),
               main = paste("Module-trait relationships - Origin"))

png("../../output/WGCNA/Treatments_Primary_heatmap2.png", 500, 1000, pointsize=20)
labeledHeatmap(Matrix = moduleTraitCor_Primary,
               xLabels = names(Traits.Origin),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = TRUE,
               colors = blueWhiteRed(50),
               textMatrix = d0.PRIMARYTreatments.matrix,
               setStdMargins = FALSE,
               cex.text = 1,
               zlim = c(-1,1),
               main = paste("Module-trait relationships - Origin"))
dev.off()
## pdf 
##   2

Module eigengene - MEs boxplots by treatment group

library(reshape2) 
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
MEs_table <- MEs # new table for plotting 
MEs_table$Coral_ID <- rownames(MEs) # call rows as coolumn to merge with treatment data
MEsPlotting <- merge(coldata, MEs_table, by = 'Coral_ID') %>% dplyr::select(-c("Colony")) # merge
#MEsPlotting <- MEsPlotting[,-c(2)] # ommit the phys data to just plot the module colors 
MEsPlotting_melt <- melt(MEsPlotting, id=c('Coral_ID', 'Origin', 'Treatment'))
#plot it

ggplot(MEsPlotting_melt, aes(x=Origin, y=value, fill = factor(Origin), shape=Origin)) +
  geom_boxplot(aes(middle = mean(value)), position=position_dodge(0.8), outlier.size = 0, alpha = 0.5) + 
  stat_summary(fun = mean, color = "black", position = position_dodge(0.75),
               geom = "point", shape = 19, size = 3,
               show.legend = FALSE) +
  ylab("ModuleEigengene") +
  ylim(-0.5,0.5) +
  scale_color_manual(values=c("#56B4E9","#D55E00")) +
  geom_hline(yintercept=0, linetype='dotted', col = 'black', size = 1)+
  theme_classic() +
  theme(legend.position = "none") +
  facet_wrap(~variable)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_summary()`).

png("../../output/WGCNA/ME_Boxplot.png", 600, 1000, pointsize=20)

ggplot(MEsPlotting_melt, aes(x=Origin, y=value, fill = factor(Origin), shape=Origin)) +
  geom_boxplot(aes(middle = mean(value)), position=position_dodge(0.8), outlier.size = 0, alpha = 0.5) + 
  stat_summary(fun = mean, color = "black", position = position_dodge(0.75),
               geom = "point", shape = 19, size = 3,
               show.legend = FALSE) +
  ylab("ModuleEigengene") +
  ylim(-0.5,0.5) +
  scale_color_manual(values=c("#56B4E9","#D55E00")) +
  geom_hline(yintercept=0, linetype='dotted', col = 'black', size = 1)+
  theme_classic() +
  theme(legend.position = "none") +
  facet_wrap(~variable)
## Warning: Removed 1 row containing non-finite outside the scale range (`stat_boxplot()`).
## Removed 1 row containing non-finite outside the scale range (`stat_summary()`).
dev.off()
## quartz_off_screen 
##                 2

Module trait correlation

Load in phys for all traits

d <- read.csv("../../data/Heron pHi coral physiology and respirometry R_reduced.csv", strip.white=T)

d$Origin <- as.factor(d$Origin)
d$Treatment <- as.factor(d$Treatment)
d$Origin <- as.numeric(d$Origin)
d$Treatment <- as.numeric(d$Treatment)
d$Colony <- as.factor(d$Colony)
d$Coral.ID <- as.factor(d$Coral.ID)
d <- plyr::rename(d, c("Coral.ID"="Coral_ID"))
d <- d[-c(7, 8),]#here we remove the two samples that are behaving weirdly
d <- subset(d, select = -c(Coral_ID))

Correlations of traits with eigengenes

moduleTraitCor_both = cor(MEs, d, use = "p")
## Warning in storage.mode(y) <- "double": NAs introduced by coercion
moduleTraitPvalue_both = corPvalueStudent(moduleTraitCor_both, nSamples)
Colors=sub("ME","",names(MEs))
moduleTraitTree = hclust(dist(t(moduleTraitCor_both)), method = "average")
plot(moduleTraitTree, main = "Clustering based on module-trait correlation", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)

### Heatmaps

#sizeGrWindow(20,20)
# Will display correlations and their p-values
both.matrix <-  paste(signif(moduleTraitCor_both, 3), "\n(",
                                       signif(moduleTraitPvalue_both, 3), ")", sep = "")
#dim(textMatrix) == dim(moduleTraitCor_treatonly)
sizeGrWindow(20,20)
labeledHeatmap(Matrix = moduleTraitCor_both,
               xLabels = names(d),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = TRUE,
               colors = blueWhiteRed(50),
               textMatrix = both.matrix,
               setStdMargins = FALSE,
               cex.text = 1,
               zlim = c(-1,1),
               main = paste("Module-trait relationships - Origin and Treatment"))

par(mar = c(10, 9.5, 5, 3));
# Display the correlation values within a heatmap plot
png("../../output/WGCNA/Both_with phys and pHi_heatmap_new.png", 2000, 2000, pointsize=20)
labeledHeatmap(Matrix = moduleTraitCor_both,
               xLabels = names(d),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = TRUE,
               colors = blueWhiteRed(50),
               textMatrix = both.matrix,
               setStdMargins = FALSE,
               cex.text = 1,
               zlim = c(-1,1),
               main = paste("Module-trait relationships - Origin and Treatment"))
dev.off()
## pdf 
##   2

Complex heatmap

library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.16.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##     genomic data. Bioinformatics 2016.
## 
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
##    in the original pheatmap() are identically supported in the new function. You 
##    can still use the original function by explicitly calling pheatmap::pheatmap().
## 
## Attaching package: 'ComplexHeatmap'
## The following object is masked from 'package:pheatmap':
## 
##     pheatmap
## The following object is masked from 'package:genefilter':
## 
##     dist2
#bold sig p-values
#dendrogram with WGCNA MEtree cut-off
#colored y-axis
#Create list of pvalues for eigengene correlation with specific life stages
heatmappval <- signif(moduleTraitPvalue_both, 1)
#Make list of heatmap row colors
htmap.colors <- rownames(moduleTraitCor_both)
htmap.colors <- gsub("ME", "", htmap.colors)

moduleTraitCor_both <- as.matrix(moduleTraitCor_both)

row_dend = as.dendrogram(hclust(dist(moduleTraitCor_both), method = "average"))
row_dend = rotate(row_dend,c(1:3,5,4,8:10, 12, 11, 14:15 ,13 ,6:7))
column_dend = as.dendrogram(hclust(dist(t(moduleTraitCor_both)), method = "average"))
column_dend = rotate(column_dend,c(1:2,12:17,7:11,3:6))

pdf("../../output/WGCNA/Both_with phys and pHi_heatmap_new_row_clust.pdf", height = 40, width = 40)
ht=Heatmap(moduleTraitCor_both, name = "Eigengene", column_title = "Module-Trait Eigengene Correlation", 
        col = blueWhiteRed(50), 
        row_names_side = "left",
        cluster_rows = row_dend,
        row_dend_width = unit(2, "in"),
        width = unit(20, "in"), height = unit(20, "in"), 
        cluster_columns = column_dend,
        column_dend_height = unit(1.25, "in"),
        row_gap = unit(2.5, "mm"), border = TRUE,
        cell_fun = function(j, i, x, y, w, h, col) {
        if(heatmappval[i, j] < 0.05) {
            grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 25, fontface = "bold"))
        }
        else {
            #grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 25, fontface = "plain"))
        }},
        column_names_gp =  gpar(fontsize = 30), row_names_gp = gpar(fontsize = 25, alpha = 0.75, border = TRUE, fill = htmap.colors))
draw(ht)
dev.off()
## quartz_off_screen 
##                 2

Module eigengene - MEs boxplots by treatment group

library(reshape2) 
MEs_table <- MEs # new table for plotting 
MEs_table$Coral_ID <- rownames(MEs) # call rows as coolumn to merge with treatment data
MEsPlotting <- merge(coldata, MEs_table, by = 'Coral_ID') %>% dplyr::select(-c("Colony")) # merge
#MEsPlotting <- MEsPlotting[,-c(2)] # ommit the phys data to just plot the module colors 
MEsPlotting_melt <- melt(MEsPlotting, id=c('Coral_ID', 'Origin', 'Treatment'))
#plot it
png("../../output/WGCNA/ME_Boxplot_Treatment.png", 600, 1000, pointsize=20)
ME_Boxplot_Treatment = ggplot(MEsPlotting_melt, aes(x=Treatment, y=value, fill = factor(Treatment), shape=Treatment)) +
  geom_boxplot(aes(middle = mean(value)), position=position_dodge(0.8), outlier.size = 0, alpha = 0.5) + 
  stat_summary(fun.y = mean, color = "black", position = position_dodge(0.75),
               geom = "point", shape = 19, size = 3,
               show.legend = FALSE) +
  ylab("ModuleEigengene") +
  ylim(-0.5,0.5) +
  scale_color_manual(values=c("#56B4E9","#D55E00")) +
  geom_hline(yintercept=0, linetype='dotted', col = 'black', size = 1)+
  theme_classic() +
  theme(legend.position = "none") +
  facet_wrap(~variable)
## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ME_Boxplot_Treatment
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_summary()`).
dev.off()
## quartz_off_screen 
##                 2
ME_Boxplot_Treatment
## Warning: Removed 1 row containing non-finite outside the scale range (`stat_boxplot()`).
## Removed 1 row containing non-finite outside the scale range (`stat_summary()`).

Plots based off of Strader & Quigley 2022: https://www.nature.com/articles/s41467-022-32217-z

# View module eigengene data and make dataframe for Strader plots.  
head(MEs)
##            MEcyan     MEblack     MEsalmon      MEpink       MEred    MEgreen
## RF13B 0.023045712  0.03573889 -0.029662215 -0.06472556 -0.04120489 0.06229807
## RF13D 0.010859993  0.01671523  0.049450663  0.05677112  0.01828699 0.04001668
## RF14B 0.008973815 -0.14656653  0.103027496 -0.02275044  0.02758936 0.09009377
## RF14C 0.055555248 -0.14441784  0.001006515 -0.02644198  0.05978645 0.12396450
## RF15B 0.066984890  0.36136333  0.307372343  0.29775735  0.26842429 0.07146349
## RF15D 0.061369555  0.34069144  0.345022730  0.26757515  0.24009722 0.12177161
##       MEmidnightblue       MEblue    MEbrown     MEpurple MEturquoise
## RF13B   -0.024549173  0.005857733 0.08621400  0.015166907 -0.04998745
## RF13D   -0.018376306  0.044219420 0.07779582  0.004543621 -0.01913715
## RF14B   -0.022880675 -0.083031859 0.24394568  0.293751818 -0.09318565
## RF14C    0.006327375  0.035898347 0.23139138  0.250251479 -0.14770968
## RF15B    0.151329333  0.033563498 0.19519443 -0.158817259  0.03523276
## RF15D    0.259392975 -0.059028765 0.13060026 -0.131872934  0.07878675
##              MEtan    MEmagenta  MElightcyan    MEyellow
## RF13B  0.004223004  0.040533178  0.004539119  0.06392860
## RF13D -0.066235131  0.004749849 -0.048921622 -0.02799559
## RF14B  0.244085251 -0.055936540  0.007711938  0.08763697
## RF14C  0.166964022 -0.015987095 -0.069472828  0.02895714
## RF15B -0.248831159 -0.362117043 -0.344930466 -0.38286220
## RF15D -0.142520988 -0.295093407 -0.167570877 -0.14015032
names(MEs)
##  [1] "MEcyan"         "MEblack"        "MEsalmon"       "MEpink"        
##  [5] "MEred"          "MEgreen"        "MEmidnightblue" "MEblue"        
##  [9] "MEbrown"        "MEpurple"       "MEturquoise"    "MEtan"         
## [13] "MEmagenta"      "MElightcyan"    "MEyellow"
Strader_MEs <- MEs
Strader_MEs$Treatment <- coldata_outrm$Treatment
Strader_MEs$Origin <- coldata_outrm$Origin
Strader_MEs$Colony <- d$Colony
Strader_MEs$Coral_ID <- rownames(Strader_MEs)
head(Strader_MEs)
##            MEcyan     MEblack     MEsalmon      MEpink       MEred    MEgreen
## RF13B 0.023045712  0.03573889 -0.029662215 -0.06472556 -0.04120489 0.06229807
## RF13D 0.010859993  0.01671523  0.049450663  0.05677112  0.01828699 0.04001668
## RF14B 0.008973815 -0.14656653  0.103027496 -0.02275044  0.02758936 0.09009377
## RF14C 0.055555248 -0.14441784  0.001006515 -0.02644198  0.05978645 0.12396450
## RF15B 0.066984890  0.36136333  0.307372343  0.29775735  0.26842429 0.07146349
## RF15D 0.061369555  0.34069144  0.345022730  0.26757515  0.24009722 0.12177161
##       MEmidnightblue       MEblue    MEbrown     MEpurple MEturquoise
## RF13B   -0.024549173  0.005857733 0.08621400  0.015166907 -0.04998745
## RF13D   -0.018376306  0.044219420 0.07779582  0.004543621 -0.01913715
## RF14B   -0.022880675 -0.083031859 0.24394568  0.293751818 -0.09318565
## RF14C    0.006327375  0.035898347 0.23139138  0.250251479 -0.14770968
## RF15B    0.151329333  0.033563498 0.19519443 -0.158817259  0.03523276
## RF15D    0.259392975 -0.059028765 0.13060026 -0.131872934  0.07878675
##              MEtan    MEmagenta  MElightcyan    MEyellow Treatment Origin
## RF13B  0.004223004  0.040533178  0.004539119  0.06392860  Variable   Flat
## RF13D -0.066235131  0.004749849 -0.048921622 -0.02799559    Stable   Flat
## RF14B  0.244085251 -0.055936540  0.007711938  0.08763697  Variable   Flat
## RF14C  0.166964022 -0.015987095 -0.069472828  0.02895714    Stable   Flat
## RF15B -0.248831159 -0.362117043 -0.344930466 -0.38286220  Variable   Flat
## RF15D -0.142520988 -0.295093407 -0.167570877 -0.14015032    Stable   Flat
##       Colony Coral_ID
## RF13B     13    RF13B
## RF13D     13    RF13D
## RF14B     14    RF14B
## RF14C     14    RF14C
## RF15B     15    RF15B
## RF15D     15    RF15D
Strader_MEs <- Strader_MEs%>%
  droplevels() #drop unused level
plot_MEs <- Strader_MEs%>%
  gather(., key="Module", value="Mean", 1:(ncol(Strader_MEs)-4))
head(plot_MEs)
##   Treatment Origin Colony Coral_ID Module        Mean
## 1  Variable   Flat     13    RF13B MEcyan 0.023045712
## 2    Stable   Flat     13    RF13D MEcyan 0.010859993
## 3  Variable   Flat     14    RF14B MEcyan 0.008973815
## 4    Stable   Flat     14    RF14C MEcyan 0.055555248
## 5  Variable   Flat     15    RF15B MEcyan 0.066984890
## 6    Stable   Flat     15    RF15D MEcyan 0.061369555
plot_MEs_mean <- summarySE(plot_MEs, measurevar='Mean', groupvars=c('Treatment','Origin', "Module"), na.rm=TRUE, conf.interval = 0.95)
plot_MEs_mean
##    Treatment Origin         Module  N         Mean         sd          se
## 1     Stable  Slope        MEblack 12 -0.039754264 0.06202534 0.017905174
## 2     Stable  Slope         MEblue 12 -0.110056761 0.20407118 0.058910276
## 3     Stable  Slope        MEbrown 12 -0.122307021 0.06850331 0.019775203
## 4     Stable  Slope         MEcyan 12  0.013369322 0.04588326 0.013245356
## 5     Stable  Slope        MEgreen 12 -0.057759479 0.18209748 0.052567016
## 6     Stable  Slope    MElightcyan 12  0.008978748 0.14457269 0.041734541
## 7     Stable  Slope      MEmagenta 12  0.094525394 0.07778958 0.022455917
## 8     Stable  Slope MEmidnightblue 12  0.022929424 0.17937023 0.051779724
## 9     Stable  Slope         MEpink 12 -0.024042545 0.13090715 0.037789639
## 10    Stable  Slope       MEpurple 12 -0.087106274 0.05005723 0.014450276
## 11    Stable  Slope          MEred 12 -0.084292611 0.12967533 0.037434044
## 12    Stable  Slope       MEsalmon 12  0.010043182 0.12510710 0.036115309
## 13    Stable  Slope          MEtan 12 -0.001827060 0.13790126 0.039808664
## 14    Stable  Slope    MEturquoise 12  0.137324385 0.15076081 0.043520897
## 15    Stable  Slope       MEyellow 12 -0.071252071 0.17217664 0.049703116
## 16    Stable   Flat        MEblack 11  0.028212561 0.20306332 0.061225894
## 17    Stable   Flat         MEblue 11  0.051722038 0.09691348 0.029220513
## 18    Stable   Flat        MEbrown 11  0.140572527 0.08941799 0.026960539
## 19    Stable   Flat         MEcyan 11  0.030429746 0.02490805 0.007510061
## 20    Stable   Flat        MEgreen 11  0.092491604 0.09756417 0.029416704
## 21    Stable   Flat    MElightcyan 11 -0.063292592 0.11152780 0.033626897
## 22    Stable   Flat      MEmagenta 11 -0.101329792 0.12740965 0.038415455
## 23    Stable   Flat MEmidnightblue 11  0.051785345 0.11968160 0.036085359
## 24    Stable   Flat         MEpink 11  0.068162754 0.14864751 0.044818912
## 25    Stable   Flat       MEpurple 11  0.086921017 0.18455066 0.055644117
## 26    Stable   Flat          MEred 11  0.123264555 0.12782853 0.038541753
## 27    Stable   Flat       MEsalmon 11  0.033156651 0.18594831 0.056065525
## 28    Stable   Flat          MEtan 11  0.025414739 0.16774005 0.050575528
## 29    Stable   Flat    MEturquoise 11 -0.094356803 0.11144867 0.033603039
## 30    Stable   Flat       MEyellow 11  0.018812643 0.12508447 0.037714387
## 31  Variable  Slope        MEblack 12 -0.021871855 0.05108377 0.014746614
## 32  Variable  Slope         MEblue 12  0.023848447 0.10688538 0.030855151
## 33  Variable  Slope        MEbrown 12 -0.125277907 0.04818037 0.013908476
## 34  Variable  Slope         MEcyan 12 -0.067307852 0.28331241 0.081785248
## 35  Variable  Slope        MEgreen 12 -0.076561117 0.15351021 0.044314581
## 36  Variable  Slope    MElightcyan 12  0.097493355 0.16605295 0.047935357
## 37  Variable  Slope      MEmagenta 12  0.098641298 0.10807566 0.031198757
## 38  Variable  Slope MEmidnightblue 12 -0.090179594 0.16882156 0.048734588
## 39  Variable  Slope         MEpink 12 -0.065797094 0.17240455 0.049768905
## 40  Variable  Slope       MEpurple 12 -0.068504434 0.03317717 0.009577425
## 41  Variable  Slope          MEred 12 -0.099034919 0.10326861 0.029811079
## 42  Variable  Slope       MEsalmon 12 -0.087344334 0.13762115 0.039727804
## 43  Variable  Slope          MEtan 12 -0.027294229 0.13771990 0.039756309
## 44  Variable  Slope    MEturquoise 12  0.024446994 0.09902844 0.028587048
## 45  Variable  Slope       MEyellow 12  0.035223379 0.14284911 0.041236987
## 46  Variable   Flat        MEblack 11  0.039015933 0.21596847 0.065116943
## 47  Variable   Flat         MEblue 11  0.042323395 0.11000329 0.033167240
## 48  Variable   Flat        MEbrown 11  0.129520122 0.08552768 0.025787566
## 49  Variable   Flat         MEcyan 11  0.028412288 0.03483993 0.010504633
## 50  Variable   Flat        MEgreen 11  0.054039955 0.06651343 0.020054555
## 51  Variable   Flat    MElightcyan 11 -0.052858793 0.12351732 0.037241874
## 52  Variable   Flat      MEmagenta 11 -0.109397509 0.13192715 0.039777533
## 53  Variable   Flat MEmidnightblue 11  0.021578478 0.07129269 0.021495555
## 54  Variable   Flat         MEpink 11  0.029844125 0.11991092 0.036154503
## 55  Variable   Flat       MEpurple 11  0.082836119 0.17714043 0.053409849
## 56  Variable   Flat          MEred 11  0.076729115 0.09973119 0.030070085
## 57  Variable   Flat       MEsalmon 11  0.051171879 0.11835171 0.035684383
## 58  Variable   Flat          MEtan 11  0.006353940 0.16860756 0.050837094
## 59  Variable   Flat    MEturquoise 11 -0.082121065 0.10597943 0.031954002
## 60  Variable   Flat       MEyellow 11  0.020491385 0.14445032 0.043553411
##            ci
## 1  0.03940902
## 2  0.12966064
## 3  0.04352493
## 4  0.02915283
## 5  0.11569922
## 6  0.09185711
## 7  0.04942514
## 8  0.11396640
## 9  0.08317443
## 10 0.03180484
## 11 0.08239177
## 12 0.07948926
## 13 0.08761828
## 14 0.09578885
## 15 0.10939582
## 16 0.13641979
## 17 0.06510736
## 18 0.06007182
## 19 0.01673346
## 20 0.06554450
## 21 0.07492540
## 22 0.08559497
## 23 0.08040319
## 24 0.09986276
## 25 0.12398282
## 26 0.08587638
## 27 0.12492177
## 28 0.11268930
## 29 0.07487224
## 30 0.08403289
## 31 0.03245708
## 32 0.06791173
## 33 0.03061235
## 34 0.18000812
## 35 0.09753574
## 36 0.10550501
## 37 0.06866800
## 38 0.10726410
## 39 0.10954062
## 40 0.02107977
## 41 0.06561374
## 42 0.08744031
## 43 0.08750305
## 44 0.06291967
## 45 0.09076200
## 46 0.14508959
## 47 0.07390122
## 48 0.05745828
## 49 0.02340578
## 50 0.04468433
## 51 0.08298007
## 52 0.08862987
## 53 0.04789508
## 54 0.08055725
## 55 0.11900456
## 56 0.06700033
## 57 0.07950976
## 58 0.11327210
## 59 0.07119795
## 60 0.09704305

By Origin, Plot mean module eigengene for each module.

png("../../output/WGCNA/plot_MEs_fig_Origin.png", 600, 1000, pointsize=20)
plot_MEs_fig_Origin <- ggplot(plot_MEs, aes(y=Mean, x=Origin, color=Origin, fill=Origin))+ 
  geom_point(alpha=0.8,position=position_jitterdodge(0.05))+
  geom_boxplot(alpha=0.4, color="black", outlier.shape = NA)+
  #geom_smooth(method=gam)+
  geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
  facet_wrap(~Module)+
  #scale_color_brewer(palette = "RdBu", direction=-1)+
  #scale_fill_brewer(palette = "RdBu", direction=-1)+
  #scale_x_continuous(expression(Temperature~(degree~C)), limits=c(23,38),breaks=c(24,26,28,30,32,34,36,38),expand = c(0.005, 0.005))+
  scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.5, 0.5), expand = c(0.005, 0.005))+  
  theme_classic()+
  theme(axis.text.x=element_text(vjust=0.5, size=12),
        axis.text.y=element_text(vjust=0.5, size=12),
        axis.title.x=element_text(size=12),
        axis.title.y=element_text(size=12),
        legend.text = element_text(vjust=0.5, size=12),
        panel.background= element_rect(fill=NA, color='black'),
        #legend.position= c(0.89,0.89),
        strip.text = element_text(vjust=0.5, size=12))
plot_MEs_fig_Origin
dev.off()
## quartz_off_screen 
##                 2

By Treatment, Plot mean module eigengene for each module.

png("../../output/WGCNA/plot_MEs_fig_Treatment.png", 600, 1000, pointsize=20)

plot_MEs_fig_Treatment <- ggplot(plot_MEs, aes(y=Mean, x=Treatment, color=Treatment, fill=Treatment))+ 
  geom_point(alpha=0.8,position=position_jitterdodge(0.05))+
  geom_boxplot(alpha=0.4, color="black", outlier.shape = NA)+
  #geom_smooth(method=gam)+
  geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
  facet_wrap(~Module)+
  scale_color_manual("Treatment", values=c("Stable"= 'lightcyan',"Variable"='orange'))+
  scale_fill_manual("Treatment", values=c("Stable"= 'lightcyan',"Variable"='orange'))+
  #scale_x_continuous(expression(Temperature~(degree~C)), limits=c(23,38),breaks=c(24,26,28,30,32,34,36,38),expand = c(0.005, 0.005))+
  scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.5, 0.5), expand = c(0.005, 0.005))+  
  theme_classic()+
  theme(axis.text.x=element_text(vjust=0.5, size=12),
        axis.text.y=element_text(vjust=0.5, size=12),
        axis.title.x=element_text(size=12),
        axis.title.y=element_text(size=12),
        legend.text = element_text(vjust=0.5, size=12),
        panel.background= element_rect(fill=NA, color='black'),
        #legend.position= c(0.89,0.89),
        strip.text = element_text(vjust=0.5, size=12))

plot_MEs_fig_Treatment
dev.off()
## quartz_off_screen 
##                 2

By Treatment and Origin, Plot mean module eigengene for each module.

png("../../output/WGCNA/plot_MEs_fig_Trt_Orig.png", 600, 1000, pointsize=20)

plot_MEs_fig_Trt_Orig <- ggplot(plot_MEs, aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+ 
  geom_point(alpha=0.8,position=position_jitterdodge(0.05))+
  geom_boxplot(alpha=0.4, color="black", outlier.shape = NA)+
  #geom_smooth(method=gam)+
  geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
  facet_wrap(~Module)+
  scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  #scale_x_continuous(expression(Temperature~(degree~C)), limits=c(23,38),breaks=c(24,26,28,30,32,34,36,38),expand = c(0.005, 0.005))+
  scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.5, 0.5), expand = c(0.005, 0.005))+  
  theme_classic()+
  theme(axis.text.x=element_text(vjust=0.5, size=12),
        axis.text.y=element_text(vjust=0.5, size=12),
        axis.title.x=element_text(size=12),
        axis.title.y=element_text(size=12),
        legend.text = element_text(vjust=0.5, size=12),
        panel.background= element_rect(fill=NA, color='black'),
        #legend.position= c(0.89,0.89),
        strip.text = element_text(vjust=0.5, size=12))
plot_MEs_fig_Trt_Orig
dev.off()
## quartz_off_screen 
##                 2

geneModuleMembership - geneTraitSignificance - GSPvalue

# Gene relationship to trait and important modules: 
# Gene Significance and Module Membership

# We quantify associations of individual genes with our trait of interest (TAOC)

# names (colors) of the modules
modNames = substring(names(MEs), 3) # name all the modules, from 3rd character on (first two are ME)
modNames
##  [1] "cyan"         "black"        "salmon"       "pink"         "red"         
##  [6] "green"        "midnightblue" "blue"         "brown"        "purple"      
## [11] "turquoise"    "tan"          "magenta"      "lightcyan"    "yellow"
geneModuleMembership = as.data.frame(cor(vst2, MEs, use = "p"));
head(geneModuleMembership)
##                                                MEcyan     MEblack    MEsalmon
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1  0.10438549 -0.12235265 -0.47423904
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 -0.16310820  0.21739302 -0.14153773
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a  0.08339074 -0.04725328 -0.33667736
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b  0.09142923  0.15329361  0.10385419
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1   0.12842828 -0.10929165  0.07806305
## Pocillopora_acuta_HIv2___RNAseq.g682.t1   -0.02392754 -0.03296640  0.17552549
##                                               MEpink       MEred     MEgreen
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 -0.3584060 -0.06200274  0.20465795
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 -0.2594287 -0.31902027 -0.01684585
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a  0.1186510  0.39184956  0.44854568
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b  0.5825581  0.82234053  0.35222463
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1   0.2478598  0.35321199  0.60087070
## Pocillopora_acuta_HIv2___RNAseq.g682.t1   -0.4538716 -0.45358620 -0.31698652
##                                           MEmidnightblue      MEblue    MEbrown
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1    -0.21979419  0.49696801  0.2454424
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1     0.08185613  0.08994886 -0.3855597
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a     0.14471353  0.48376147  0.3927873
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b     0.05633416  0.51529570  0.6432814
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1      0.44897525  0.18646978  0.4855880
## Pocillopora_acuta_HIv2___RNAseq.g682.t1      -0.10512744 -0.50077743 -0.2318526
##                                              MEpurple MEturquoise       MEtan
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1  0.35441847  -0.5622244  0.15275154
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 -0.41677235   0.1495377 -0.34664520
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a  0.34074711  -0.5934280  0.02119366
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b  0.33142341  -0.5855358 -0.20954032
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1   0.39122310  -0.3553428  0.10450615
## Pocillopora_acuta_HIv2___RNAseq.g682.t1   -0.04828805   0.4289215  0.41719038
##                                             MEmagenta MElightcyan    MEyellow
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1  0.07322422  0.28404105  0.53716214
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1  0.16236103  0.15377483  0.01771476
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a -0.21754628 -0.06376141  0.33701433
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b -0.66545958 -0.41176150  0.08707934
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1  -0.38686031 -0.34369924 -0.01298406
## Pocillopora_acuta_HIv2___RNAseq.g682.t1    0.26807670  0.31248662  0.01298851
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");


# Flat treatment group
Flat = as.data.frame(Traits.Origin$Flat); # Define variable containing the desired column 
names(Flat) = "Flat"
Flat_geneTraitSignificance = as.data.frame(cor(vst2, Flat, use = "p"));
Flat_GSPvalue = as.data.frame(corPvalueStudent(as.matrix(Flat_geneTraitSignificance), nSamples));
names(Flat_geneTraitSignificance) = paste("GS.", names(Flat), sep=""); # MA_geneTraitSignificance - pearsons correlation between reads and the MA grop
names(Flat_GSPvalue) = paste("p.GS.", names(Flat), sep=""); # corPvalueStudent

# Slope treatment group
Slope = as.data.frame(Traits.Origin$Slope); # Define variable containing the desired column 
names(Slope) = "Slope"
Slope_geneTraitSignificance = as.data.frame(cor(vst2, Slope, use = "p"));
Slope_GSPvalue = as.data.frame(corPvalueStudent(as.matrix(Slope_geneTraitSignificance), nSamples));
names(Slope_geneTraitSignificance) = paste("GS.", names(Slope), sep=""); # MA_geneTraitSignificance - pearsons correlation between reads and the MA grop
names(Slope_GSPvalue) = paste("p.GS.", names(Slope), sep=""); # corPvalueStudent

Call annotation data to get module gene data (prep for downstream GO)

library(rtracklayer)
gff <- rtracklayer::import("../../data/Pocillopora_acuta_HIv2.genes_fixed.gff3")
gff <- as.data.frame(gff)
dim(gff) # 478988     13
## [1] 478988     13
names(gff) 
##  [1] "seqnames"      "start"         "end"           "width"        
##  [5] "strand"        "source"        "type"          "score"        
##  [9] "phase"         "ID"            "transcript_id" "gene_id"      
## [13] "Parent"
transcripts <- subset(gff, type == "transcript")
transcripts_gr <- makeGRangesFromDataFrame(transcripts, keep.extra.columns=TRUE) #extract length information
transcript_lengths <- width(transcripts_gr) #isolate length of each gene
seqnames <- transcripts_gr$ID #extract list of gene id 
lengths <- cbind(seqnames, transcript_lengths)
lengths <- as.data.frame(lengths) #convert to data frame
eggnog <- read.delim("../../data/Pocillopora_acuta_HIv2.genes.EggNog_results.txt")
eggnog <- plyr::rename(eggnog, c("X.query"="gene_id"))
gogene <- merge(transcripts, eggnog, by=c("gene_id"), all=T)
colnames(gogene)
##  [1] "gene_id"        "seqnames"       "start"          "end"           
##  [5] "width"          "strand"         "source"         "type"          
##  [9] "score.x"        "phase"          "ID"             "transcript_id" 
## [13] "Parent"         "seed_ortholog"  "evalue"         "score.y"       
## [17] "eggNOG_OGs"     "max_annot_lvl"  "COG_category"   "Description"   
## [21] "Preferred_name" "GOs"            "EC"             "KEGG_ko"       
## [25] "KEGG_Pathway"   "KEGG_Module"    "KEGG_Reaction"  "KEGG_rclass"   
## [29] "BRITE"          "KEGG_TC"        "CAZy"           "BiGG_Reaction" 
## [33] "PFAMs"
probes = rownames(cts_filtpoa)
probes2annot = match(probes, gogene$gene_id)
# The following is the number or probes without annotation:
sum(is.na(probes2annot))
## [1] 0
# Should return 0.

BUILD GENE INFO DATAFRAMES

length(probes)
## [1] 9012
length(gogene$seqnames[probes2annot])
## [1] 9012
length(gogene$seqnames[probes2annot])
## [1] 9012
# Create the starting data frame
names(Flat_geneTraitSignificance)
## [1] "GS.Flat"
names(Flat_GSPvalue)
## [1] "p.GS.Flat"
geneInfo_GROUPS = data.frame(substanceBXH = probes,
                                  geneSymbol = gogene$seqnames[probes2annot],
                                  moduleColor = moduleColors,
                                  GO.terms = gogene$GOs[probes2annot],
                                  GO.description = gogene$Description[probes2annot],
                                  Flat_geneTraitSignificance, Slope_geneTraitSignificance, # call this specific to the module and trait of interest
                                  Flat_GSPvalue,  Slope_GSPvalue)              # call this specific to the module and trait of interest
head(geneInfo_GROUPS)
##                                                                        substanceBXH
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 Pocillopora_acuta_HIv2___RNAseq.g27841.t1
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 Pocillopora_acuta_HIv2___RNAseq.g14011.t1
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a Pocillopora_acuta_HIv2___RNAseq.g7479.t3a
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b Pocillopora_acuta_HIv2___RNAseq.g7479.t3b
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1   Pocillopora_acuta_HIv2___RNAseq.g3340.t1
## Pocillopora_acuta_HIv2___RNAseq.g682.t1     Pocillopora_acuta_HIv2___RNAseq.g682.t1
##                                                                     geneSymbol
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1   Pocillopora_acuta_HIv2___Sc0000000
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1   Pocillopora_acuta_HIv2___Sc0000005
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a Pocillopora_acuta_HIv2___xfSc0000001
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b Pocillopora_acuta_HIv2___xfSc0000001
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1  Pocillopora_acuta_HIv2___xfSc0000006
## Pocillopora_acuta_HIv2___RNAseq.g682.t1   Pocillopora_acuta_HIv2___xfSc0000002
##                                           moduleColor
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1      yellow
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1   turquoise
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a        blue
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b         red
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1        green
## Pocillopora_acuta_HIv2___RNAseq.g682.t1     turquoise
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       GO.terms
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      GO:0000003,GO:0003006,GO:0003674,GO:0003824,GO:0004672,GO:0004674,GO:0004693,GO:0005575,GO:0005622,GO:0005623,GO:0005634,GO:0005737,GO:0005813,GO:0005815,GO:0005856,GO:0006464,GO:0006468,GO:0006793,GO:0006796,GO:0006807,GO:0007154,GO:0007165,GO:0007548,GO:0008150,GO:0008152,GO:0009987,GO:0015630,GO:0016301,GO:0016310,GO:0016740,GO:0016772,GO:0016773,GO:0019538,GO:0022414,GO:0023052,GO:0032502,GO:0036211,GO:0043170,GO:0043226,GO:0043227,GO:0043228,GO:0043229,GO:0043231,GO:0043232,GO:0043412,GO:0044237,GO:0044238,GO:0044260,GO:0044267,GO:0044422,GO:0044424,GO:0044430,GO:0044446,GO:0044464,GO:0050789,GO:0050794,GO:0050896,GO:0051716,GO:0051726,GO:0065007,GO:0071704,GO:0097472,GO:0140096,GO:1901564
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1                                                        GO:0000118,GO:0000123,GO:0003674,GO:0003676,GO:0003677,GO:0003712,GO:0003713,GO:0005488,GO:0005575,GO:0005622,GO:0005623,GO:0005634,GO:0005654,GO:0005730,GO:0006325,GO:0006338,GO:0006355,GO:0006357,GO:0006464,GO:0006473,GO:0006475,GO:0006807,GO:0006996,GO:0008150,GO:0008152,GO:0009889,GO:0009891,GO:0009893,GO:0009987,GO:0010468,GO:0010556,GO:0010557,GO:0010604,GO:0016043,GO:0016569,GO:0016570,GO:0016573,GO:0016604,GO:0016607,GO:0018193,GO:0018205,GO:0018393,GO:0018394,GO:0019219,GO:0019222,GO:0019538,GO:0030914,GO:0031248,GO:0031323,GO:0031325,GO:0031326,GO:0031328,GO:0031974,GO:0031981,GO:0032991,GO:0036211,GO:0043170,GO:0043226,GO:0043227,GO:0043228,GO:0043229,GO:0043231,GO:0043232,GO:0043233,GO:0043412,GO:0043543,GO:0043966,GO:0044237,GO:0044238,GO:0044260,GO:0044267,GO:0044422,GO:0044424,GO:0044428,GO:0044446,GO:0044451,GO:0044464,GO:0045935,GO:0048518,GO:0048522,GO:0050789,GO:0050794,GO:0051171,GO:0051173,GO:0051252,GO:0051254,GO:0051276,GO:0060255,GO:0065007,GO:0070013,GO:0070461,GO:0071704,GO:0071840,GO:0080090,GO:0097159,GO:0140110,GO:1901363,GO:1901564,GO:1902493,GO:1902494,GO:1902680,GO:1903506,GO:1903508,GO:1990234,GO:2000112,GO:2001141
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 <NA>
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 <NA>
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1  GO:0003674,GO:0003824,GO:0005488,GO:0005515,GO:0005575,GO:0005622,GO:0005623,GO:0005737,GO:0005739,GO:0005740,GO:0005743,GO:0005759,GO:0006629,GO:0006720,GO:0006732,GO:0006733,GO:0006743,GO:0006744,GO:0006950,GO:0007275,GO:0007399,GO:0008150,GO:0008152,GO:0008299,GO:0008610,GO:0009058,GO:0009108,GO:0009987,GO:0010941,GO:0014016,GO:0014019,GO:0016020,GO:0016043,GO:0016740,GO:0016765,GO:0019866,GO:0019898,GO:0022008,GO:0022607,GO:0030154,GO:0031090,GO:0031312,GO:0031314,GO:0031966,GO:0031967,GO:0031974,GO:0031975,GO:0032501,GO:0032502,GO:0032991,GO:0033554,GO:0042180,GO:0042181,GO:0042981,GO:0043066,GO:0043067,GO:0043069,GO:0043226,GO:0043227,GO:0043229,GO:0043231,GO:0043233,GO:0043933,GO:0044085,GO:0044237,GO:0044238,GO:0044249,GO:0044255,GO:0044281,GO:0044283,GO:0044422,GO:0044424,GO:0044425,GO:0044429,GO:0044444,GO:0044446,GO:0044455,GO:0044464,GO:0046982,GO:0046983,GO:0048468,GO:0048519,GO:0048523,GO:0048699,GO:0048731,GO:0048856,GO:0048863,GO:0048864,GO:0048869,GO:0050789,GO:0050794,GO:0050896,GO:0051186,GO:0051188,GO:0051259,GO:0051262,GO:0051290,GO:0051291,GO:0051716,GO:0060548,GO:0065003,GO:0065007,GO:0070013,GO:0071704,GO:0071840,GO:0098780,GO:1901576,GO:1901661,GO:1901663,GO:1902494,GO:1990234
## Pocillopora_acuta_HIv2___RNAseq.g682.t1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   <NA>
##                                                                                      GO.description
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 cyclin-dependent protein serine/threonine kinase activity
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1          TAF6-like RNA polymerase II, p300 CBP-associated
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a                                                      <NA>
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b                                                      <NA>
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1                  trans-hexaprenyltranstransferase activity
## Pocillopora_acuta_HIv2___RNAseq.g682.t1                                                        <NA>
##                                              GS.Flat   GS.Slope    p.GS.Flat
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1  0.2031928 -0.2031928 1.756178e-01
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 -0.3260066  0.3260066 2.703391e-02
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a  0.3491137 -0.3491137 1.740688e-02
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b  0.6161971 -0.6161971 5.140317e-06
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1   0.3860958 -0.3860958 8.044440e-03
## Pocillopora_acuta_HIv2___RNAseq.g682.t1   -0.1301845  0.1301845 3.885043e-01
##                                             p.GS.Slope
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 1.756178e-01
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 2.703391e-02
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a 1.740688e-02
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b 5.140317e-06
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1  8.044440e-03
## Pocillopora_acuta_HIv2___RNAseq.g682.t1   3.885043e-01
modOrder = order(-abs(cor(MEs, Flat, use = "p"))); # order by the strength of the correlation between module and trait values for each sample

for (mod in 1:ncol(geneModuleMembership)) # Add module membership information in the chosen order
{
  oldNames = names(geneInfo_GROUPS)
  geneInfo_GROUPS = data.frame(geneInfo_GROUPS, geneModuleMembership[, modOrder[mod]], 
                                    MMPvalue[, modOrder[mod]]);
  names(geneInfo_GROUPS) = c(oldNames, paste("A.", modNames[modOrder[mod]], sep=""),
                                  paste("p.A.", modNames[modOrder[mod]], sep=""))
}

Order the genes in the geneInfo variable first by module color, then by geneTraitSignificance

geneOrder = order(geneInfo_GROUPS$moduleColor, -abs(geneInfo_GROUPS$GS.Flat));
geneInfo_GROUPS = geneInfo_GROUPS[geneOrder, ]
head(geneInfo_GROUPS)
##                                                                        substanceBXH
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 Pocillopora_acuta_HIv2___RNAseq.g11459.t1
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 Pocillopora_acuta_HIv2___RNAseq.g27371.t1
## Pocillopora_acuta_HIv2___RNAseq.30217_t     Pocillopora_acuta_HIv2___RNAseq.30217_t
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 Pocillopora_acuta_HIv2___RNAseq.g29237.t1
## Pocillopora_acuta_HIv2___TS.g23133.t3         Pocillopora_acuta_HIv2___TS.g23133.t3
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 Pocillopora_acuta_HIv2___RNAseq.g12199.t1
##                                                                     geneSymbol
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1   Pocillopora_acuta_HIv2___Sc0000013
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 Pocillopora_acuta_HIv2___xfSc0000042
## Pocillopora_acuta_HIv2___RNAseq.30217_t   Pocillopora_acuta_HIv2___xpSc0000388
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 Pocillopora_acuta_HIv2___xfSc0000022
## Pocillopora_acuta_HIv2___TS.g23133.t3     Pocillopora_acuta_HIv2___xfSc0000008
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 Pocillopora_acuta_HIv2___xfSc0000009
##                                           moduleColor
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1       black
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1       black
## Pocillopora_acuta_HIv2___RNAseq.30217_t         black
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1       black
## Pocillopora_acuta_HIv2___TS.g23133.t3           black
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1       black
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   GO.terms
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             <NA>
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             <NA>
## Pocillopora_acuta_HIv2___RNAseq.30217_t                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  -
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             <NA>
## Pocillopora_acuta_HIv2___TS.g23133.t3     GO:0000003,GO:0001101,GO:0001553,GO:0001868,GO:0001869,GO:0002020,GO:0002376,GO:0002437,GO:0002438,GO:0002526,GO:0002576,GO:0002673,GO:0002682,GO:0002683,GO:0002697,GO:0002698,GO:0002920,GO:0002921,GO:0003006,GO:0003674,GO:0004857,GO:0004866,GO:0004867,GO:0005096,GO:0005102,GO:0005488,GO:0005515,GO:0005575,GO:0005576,GO:0005615,GO:0005622,GO:0005623,GO:0005737,GO:0005829,GO:0006807,GO:0006810,GO:0006887,GO:0006950,GO:0006952,GO:0006953,GO:0006954,GO:0006955,GO:0007275,GO:0007548,GO:0007565,GO:0007566,GO:0007584,GO:0007596,GO:0007597,GO:0007599,GO:0008047,GO:0008150,GO:0008152,GO:0008406,GO:0008585,GO:0009605,GO:0009611,GO:0009719,GO:0009725,GO:0009790,GO:0009892,GO:0009966,GO:0009987,GO:0009991,GO:0010033,GO:0010037,GO:0010466,GO:0010468,GO:0010605,GO:0010629,GO:0010646,GO:0010951,GO:0010955,GO:0012505,GO:0014070,GO:0016043,GO:0016192,GO:0019222,GO:0019538,GO:0019838,GO:0019899,GO:0019955,GO:0019956,GO:0019958,GO:0019959,GO:0019966,GO:0022411,GO:0022414,GO:0022602,GO:0022607,GO:0022617,GO:0023051,GO:0030141,GO:0030154,GO:0030162,GO:0030198,GO:0030234,GO:0030414,GO:0030449,GO:0030695,GO:0031091,GO:0031093,GO:0031323,GO:0031324,GO:0031347,GO:0031348,GO:0031410,GO:0031667,GO:0031960,GO:0031974,GO:0031982,GO:0031983,GO:0032101,GO:0032268,GO:0032269,GO:0032501,GO:0032502,GO:0032504,GO:0032940,GO:0033993,GO:0034694,GO:0034695,GO:0034774,GO:0042060,GO:0042221,GO:0042493,GO:0042698,GO:0042802,GO:0042803,GO:0043062,GO:0043085,GO:0043086,GO:0043087,GO:0043120,GO:0043121,GO:0043170,GO:0043226,GO:0043227,GO:0043229,GO:0043233,GO:0043547,GO:0043933,GO:0044085,GO:0044092,GO:0044093,GO:0044238,GO:0044421,GO:0044422,GO:0044424,GO:0044433,GO:0044444,GO:0044446,GO:0044464,GO:0044703,GO:0044706,GO:0044877,GO:0045055,GO:0045088,GO:0045137,GO:0045824,GO:0045861,GO:0045916,GO:0046545,GO:0046660,GO:0046903,GO:0046983,GO:0048306,GO:0048403,GO:0048406,GO:0048511,GO:0048513,GO:0048519,GO:0048523,GO:0048545,GO:0048568,GO:0048583,GO:0048585,GO:0048608,GO:0048609,GO:0048731,GO:0048856,GO:0048863,GO:0048869,GO:0050727,GO:0050776,GO:0050777,GO:0050789,GO:0050790,GO:0050794,GO:0050817,GO:0050878,GO:0050896,GO:0051056,GO:0051171,GO:0051172,GO:0051179,GO:0051234,GO:0051246,GO:0051248,GO:0051259,GO:0051260,GO:0051262,GO:0051289,GO:0051336,GO:0051345,GO:0051346,GO:0051384,GO:0051704,GO:0052547,GO:0052548,GO:0060205,GO:0060255,GO:0060589,GO:0061134,GO:0061135,GO:0061458,GO:0065003,GO:0065007,GO:0065008,GO:0065009,GO:0070013,GO:0070542,GO:0070613,GO:0071704,GO:0071840,GO:0072347,GO:0072376,GO:0072378,GO:0080090,GO:0080134,GO:0097305,GO:0097708,GO:0098772,GO:0099503,GO:1901564,GO:1901654,GO:1901700,GO:1902531,GO:1903317,GO:1903318,GO:1990402,GO:2000257,GO:2000258
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             <NA>
##                                                             GO.description
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1                             <NA>
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1                             <NA>
## Pocillopora_acuta_HIv2___RNAseq.30217_t          GTPase activator activity
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1                             <NA>
## Pocillopora_acuta_HIv2___TS.g23133.t3     endopeptidase inhibitor activity
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1                             <NA>
##                                             GS.Flat   GS.Slope    p.GS.Flat
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.5677405 -0.5677405 3.872268e-05
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.5452501 -0.5452501 8.913481e-05
## Pocillopora_acuta_HIv2___RNAseq.30217_t   0.5127223 -0.5127223 2.693006e-04
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 0.4741063 -0.4741063 8.732422e-04
## Pocillopora_acuta_HIv2___TS.g23133.t3     0.4645548 -0.4645548 1.144263e-03
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 0.4551104 -0.4551104 1.483628e-03
##                                             p.GS.Slope   A.brown  p.A.brown
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 3.872268e-05 0.3735665 0.01055044
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 8.913481e-05 0.3655440 0.01248495
## Pocillopora_acuta_HIv2___RNAseq.30217_t   2.693006e-04 0.2255392 0.13179089
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 8.732422e-04 0.3931221 0.00687840
## Pocillopora_acuta_HIv2___TS.g23133.t3     1.144263e-03 0.3318744 0.02424532
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 1.483628e-03 0.2863948 0.05365629
##                                            A.magenta  p.A.magenta     A.red
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 -0.4390326 2.270613e-03 0.3067458
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 -0.5153303 2.474598e-04 0.2811921
## Pocillopora_acuta_HIv2___RNAseq.30217_t   -0.4587281 1.344296e-03 0.4699351
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 -0.6079115 7.434349e-06 0.5714810
## Pocillopora_acuta_HIv2___TS.g23133.t3     -0.5510449 7.231490e-05 0.3696800
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 -0.5088948 3.045180e-04 0.4181299
##                                                p.A.red A.turquoise
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 3.812803e-02 -0.10583159
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 5.835180e-02 -0.21881675
## Pocillopora_acuta_HIv2___RNAseq.30217_t   9.835954e-04 -0.09531465
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 3.350908e-05 -0.04948723
## Pocillopora_acuta_HIv2___TS.g23133.t3     1.145282e-02 -0.23930488
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 3.832660e-03 -0.03439068
##                                           p.A.turquoise     A.purple p.A.purple
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1     0.4839269 -0.001344154  0.9929264
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1     0.1440109  0.097484918  0.5192403
## Pocillopora_acuta_HIv2___RNAseq.30217_t       0.5286280 -0.242928037  0.1037909
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1     0.7439690 -0.085020600  0.5742610
## Pocillopora_acuta_HIv2___TS.g23133.t3         0.1092114 -0.051265727  0.7350959
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1     0.8205039 -0.060136362  0.6913677
##                                               A.green    p.A.green A.lightcyan
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1  0.23766284 1.117381e-01  -0.1952640
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1  0.26820482 7.151857e-02  -0.1217634
## Pocillopora_acuta_HIv2___RNAseq.30217_t    0.16609603 2.699402e-01  -0.2728921
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 -0.04119175 7.857737e-01  -0.3495440
## Pocillopora_acuta_HIv2___TS.g23133.t3      0.55430775 6.417267e-05  -0.3123909
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1  0.20017449 1.822589e-01  -0.4310812
##                                           p.A.lightcyan     A.pink     p.A.pink
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1   0.193444104 0.18853025 0.2095613030
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1   0.420174721 0.08986363 0.5525688928
## Pocillopora_acuta_HIv2___RNAseq.30217_t     0.066519359 0.37522665 0.0101840079
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1   0.017259628 0.53290057 0.0001374455
## Pocillopora_acuta_HIv2___TS.g23133.t3       0.034545563 0.25217327 0.0908943212
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1   0.002781693 0.44616937 0.0018845316
##                                                A.blue  p.A.blue  A.salmon
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1  0.03781604 0.8029651 0.4435399
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1  0.16078732 0.2857617 0.3429064
## Pocillopora_acuta_HIv2___RNAseq.30217_t    0.18424167 0.2203025 0.3152706
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1  0.07874465 0.6029359 0.4624926
## Pocillopora_acuta_HIv2___TS.g23133.t3      0.27193253 0.0675195 0.2510276
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 -0.03636981 0.8103591 0.4161561
##                                            p.A.salmon A.midnightblue
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.002019420     0.34284439
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.019652076     0.22524824
## Pocillopora_acuta_HIv2___RNAseq.30217_t   0.032827874     0.24143271
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 0.001211796     0.01104787
## Pocillopora_acuta_HIv2___TS.g23133.t3     0.092421359     0.55037413
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 0.004020175     0.33515410
##                                           p.A.midnightblue   A.black
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1     1.967570e-02 0.5029437
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1     1.323031e-01 0.4355470
## Pocillopora_acuta_HIv2___RNAseq.30217_t       1.060024e-01 0.6425922
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1     9.419093e-01 0.4207111
## Pocillopora_acuta_HIv2___TS.g23133.t3         7.410146e-05 0.4375792
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1     2.279446e-02 0.4473112
##                                              p.A.black       A.cyan   p.A.cyan
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 3.675918e-04  0.200348223 0.18187183
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 2.483421e-03  0.273605353 0.06578363
## Pocillopora_acuta_HIv2___RNAseq.30217_t   1.474083e-06 -0.009882091 0.94803020
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 3.599055e-03 -0.045947451 0.76172625
## Pocillopora_acuta_HIv2___TS.g23133.t3     2.357306e-03  0.179255763 0.23326230
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 1.828487e-03  0.149449675 0.32153753
##                                              A.yellow p.A.yellow      A.tan
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 -0.04482242 0.76739590 -0.0465158
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1  0.13687184 0.36438494 -0.1011951
## Pocillopora_acuta_HIv2___RNAseq.30217_t   -0.10080009 0.50506242 -0.4158460
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 -0.21492986 0.15144861 -0.3403418
## Pocillopora_acuta_HIv2___TS.g23133.t3      0.03639833 0.81021311 -0.3617676
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 -0.27200931 0.06743904 -0.2045885
##                                               p.A.tan
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.758866700
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.503386444
## Pocillopora_acuta_HIv2___RNAseq.30217_t   0.004050348
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 0.020648495
## Pocillopora_acuta_HIv2___TS.g23133.t3     0.013495757
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 0.172606503
write.csv(geneInfo_GROUPS, file = "../../output/WGCNA/WGCNA_ModuleMembership.csv")

Individual Module Plots again, with n = number of genes in module

nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="brown")) #942
## [1] 942
nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="red")) #425
## [1] 425
nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="turquoise")) #2558
## [1] 2558
nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="magenta")) #219
## [1] 219
Brown only - By Treatment and Origin, Plot mean module eigengene
n_mod <- nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="brown")) #942

pdf("../../output/WGCNA/plot_MEs_fig_Trt_Orig_BROWN.pdf")

plot_MEs_fig_Trt_Orig <- ggplot(filter(plot_MEs, Module =="MEbrown"), aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+ 
  geom_boxplot(alpha=0.4, color="black", outlier.shape = NA, size =1)+
  geom_point(alpha=0.8,position=position_jitterdodge(0.05),size=2.5)+
  geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
  scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.4, 0.4), expand = c(0.005, 0.005))+  
  theme_classic()+
  labs(title = "Brown Module") +
  annotate("text", x = Inf, y = Inf, label = paste("n genes =", n_mod), hjust = 1.1, vjust = 2.1, size = 6, color = "black") + 
  theme(axis.text.x=element_text(vjust=0.5, size=20),
        axis.text.y=element_text(vjust=0.5, size=20),
        axis.title.x=element_blank(),
        axis.title.y=element_text(size=20),
        legend.text = element_text(vjust=0.5, size=20),
        panel.background= element_rect(fill=NA, color='black'),
        strip.text = element_text(vjust=0.5, size=20),
        legend.position = "none") 

plot_MEs_fig_Trt_Orig
dev.off()
## quartz_off_screen 
##                 2
plot_MEs_fig_Trt_Orig

Red only - By Treatment and Origin, Plot mean module eigengene
n_mod <- nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="red")) #425

pdf("../../output/WGCNA/plot_MEs_fig_Trt_Orig_RED.pdf")

plot_MEs_fig_Trt_Orig <- ggplot(filter(plot_MEs, Module =="MEred"), aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+ 
  geom_boxplot(alpha=0.4, color="black", outlier.shape = NA, size =1)+
  geom_point(alpha=0.8,position=position_jitterdodge(0.05),size=2.5)+
  geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
  scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.4, 0.4), expand = c(0.005, 0.005))+  
  theme_classic()+
  labs(title = "Red Module")+
  annotate("text", x = Inf, y = Inf, label = paste("n genes =", n_mod), hjust = 1.1, vjust = 2.1, size = 6, color = "black") + 
  theme(axis.text.x=element_text(vjust=0.5, size=20),
        axis.text.y=element_text(vjust=0.5, size=20),
        axis.title.x=element_blank(),
        axis.title.y=element_text(size=20),
        legend.text = element_text(vjust=0.5, size=20),
        panel.background= element_rect(fill=NA, color='black'),
        strip.text = element_text(vjust=0.5, size=20),
        legend.position = "none")
plot_MEs_fig_Trt_Orig
dev.off()
## quartz_off_screen 
##                 2
plot_MEs_fig_Trt_Orig

Turquoise only - By Treatment and Origin, Plot mean module eigengene
n_mod <- nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="turquoise")) #2558

pdf("../../output/WGCNA/plot_MEs_fig_Trt_Orig_TURQUOISE.pdf")

plot_MEs_fig_Trt_Orig <- ggplot(filter(plot_MEs, Module =="MEturquoise"), aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+ 
  geom_boxplot(alpha=0.4, color="black", outlier.shape = NA, size =1)+
  geom_point(alpha=0.8,position=position_jitterdodge(0.05),size=2.5)+
  geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
  scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.4, 0.4), expand = c(0.005, 0.005))+  
  theme_classic()+
  labs(title = "Turquoise Module")+
  annotate("text", x = Inf, y = Inf, label = paste("n genes =", n_mod), hjust = 1.1, vjust = 2.1, size = 6, color = "black") + 
  theme(axis.text.x=element_text(vjust=0.5, size=20),
        axis.text.y=element_text(vjust=0.5, size=20),
        axis.title.x=element_blank(),
        axis.title.y=element_text(size=20),
        legend.text = element_text(vjust=0.5, size=20),
        panel.background= element_rect(fill=NA, color='black'),
        strip.text = element_text(vjust=0.5, size=20),
        legend.position = "none")
plot_MEs_fig_Trt_Orig
dev.off()
## quartz_off_screen 
##                 2
plot_MEs_fig_Trt_Orig

Magenta only - By Treatment and Origin, Plot mean module eigengene
n_mod <- nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="magenta")) #219

pdf("../../output/WGCNA/plot_MEs_fig_Trt_Orig_MAGENTA.pdf")

plot_MEs_fig_Trt_Orig <- ggplot(filter(plot_MEs, Module =="MEmagenta"), aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+ 
  geom_boxplot(alpha=0.4, color="black", outlier.shape = NA, size =1)+
  geom_point(alpha=0.8,position=position_jitterdodge(0.05),size=2.5)+
  geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
  scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
  scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.4, 0.4), expand = c(0.005, 0.005))+  
  theme_classic()+
  labs(title = "Magenta Module")+
  annotate("text", x = Inf, y = Inf, label = paste("n genes =", n_mod), hjust = 1.1, vjust = 2.1, size = 6, color = "black") + 
  theme(axis.text.x=element_text(vjust=0.5, size=20),
        axis.text.y=element_text(vjust=0.5, size=20),
        axis.title.x=element_blank(),
        axis.title.y=element_text(size=20),
        legend.text = element_text(vjust=0.5, size=20),
        panel.background= element_rect(fill=NA, color='black'),
        strip.text = element_text(vjust=0.5, size=20),
        legend.position = "none")
plot_MEs_fig_Trt_Orig
dev.off()
## quartz_off_screen 
##                 2
plot_MEs_fig_Trt_Orig